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"■  An  inverse  formulation  is  developed  for  solving  three-dimensional  potential 
fluid  flows  which  considers  the  magnitudes  of  the  cartesian  coordinates  x,  y,  and  z 
as  the  dependent  variables  in  the  space  defined  by  (i.e.  the  independent  variables) 
the  potential  functior  and  two  mutually  orthogonal  stream  surface  functions  whose 
intersection  defines  the  physical  space  streamlines.  This  formulation  reverses  the 
usual  role  of  the  variables.  In  this  inverse  space  irregular  boundaries,  witfr  unknown 
position  in  the  physical  space,  such  as  free  surfaces  become  plane  boundaries,  and 
the  space  of  most  potential  flow  problems  is  a  parallelepiped. 

The  basic  partial  differential  equations  resulting  from  this  formulation 
are  nonlinear  and  three  in  number.  Finite  difference  methods  are  developed  for 
solving  the  space  boundary  value  problems  simultaneously,  which  are  associated 
with  these  three  equations.  The  applicability  of  the  inverse  formulation  and  the 
numerical  solution  is  demonstrafed'hy'  obtaining  a  solution  to  the  three-dimensional, 
free  surface  flow  past  a  vertical  strut  which  extends  through  the  fluid  surface 
and  is  placed  between  channel  walls.  | 
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ABSTRACT 


An  inverse  formulation  is  developed  for  solving  three-dimensional  potential  fluid  flows 
which  considers  the  magnitudes  of  the  cartesian  coordinates  x,  y,  and  z  as  the  dependent  variables 
in  the  space  defined  by  the  potential  function  and  two  mutually  orthogonal  stream  surface 
functions  whose  intersection  defines  the  physical  space  streamlines.  This  formulation  reverses  the 
usual  role  of  the  variables.  In  this  inverse  space  irregular'  boundaries,  with  unknown  position  in  the 
physical  space,  suck  as  free  surfaces  become  plane  boundaries,  and  the  space  of  most  potential 
flow  problems  is  a  parallelepiped. 

The  basic  partial  differential  equations  resulting  from  this  formulation  are  nonlinear  and 
three  in  number.  Finite  difference  methods  are  developed  for  solving  the  space  boundary  value 
problems  simultaneously,  which  are  associated  with  these  three  equations.  The  applicability  of  the 
inverse  formulation  and  the  numerical  solution  is  demonstrated  by  obtaining  a  solution  to  the 
three-dimensional,  free  surface  flow  past  a  vertical  strut  which  extends  through  the  fluid  surface 
and  is  placed  between  channel  walls. 
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NOTATION 


t  f 

I 


any  vector  quantity 
3x/3«j> 

acceleration  of  gravity 
any  vector  quantity 
3y/  3  4> 


Jacobian  determinant 


inverse  Jacobian  determinant 

subscript  denoting  increment  in  i|>  direction 

unit  vector  in  y-direction 

subscript  denoting  increment  in  'J'*  direction 


unit  vector  in  z-direction 


D  =  depth  of  upstream  flow 

D  =  derivative  determinant 

d  =  3y/  3i{i 

e  =  3z/3i^ 

F  =  denotes  function  of 

f  =  denotes  function  of 

f  =  3y/3'Ji* 

f  =  vector  array  of  values 
G  =  denotes  function  of 

e  =  denotes  function  of 


depth  of  flow  plus  velocity  head 
denotes  function  of 
denotes  function  of 


number  of  <P  grid  planes 


number  of  tji  grid  planes 


number  of  grid  planes 


0  *  plane  coincident  with  strut 
superscript  denoting  iteration  number 


subscript  denoting  increment  in  $  direction 
unit  vector  in  x-direction 


flow  rate 


superscript  denoting  iteration  number 

vector  array  of  values 

vector  array  of  values 

vector  array  of  values 

velocity  component  in  x-direction 

velocity  component  in  y-direction 

velocity  component  in  z-direction 


over-relaxation  factor 


width  of  channel 

0 

* 

angle 

cartesian  coordinate  and  also  magnitude  thereof 

IT 

= 

3.1417 

cartesian  coordinate  and  also  magnitude  thereof 

$ 

= 

potential  function 

cartesian  coordinate  and  also  magnitude  thereof 

4> 

= 

dimensionless  ttotential  function 

direct'  •*  cosine 

Y 

= 

stream  surface  function 

direction  cosine 

V 

* 

'dimensionless  stream  surface  function 

i  -1  .  ri. 

1  3x  dy  dr 

n* 

* 

stream  surface  function 

direction  cosine 

s 

dimensionless  stream  surface  function 

r*7**ffW 


INTRODUCTION 


For  many  practical  fluid  flow  problems  in  which 
viscous  forces  are  of  minor  importance,  because  they  are 
confined  to  relatively  small  regions  of  the  flow,  inviscid 
fluid  flow  theory  yields  results  which  are  adequate  for 
most  applications.  Consequently  a  vast  amount  of  litera¬ 
ture  deals  with-  inviscid  fluid  flow  theory.  Despite  all  of 
the  effort  represented  by  this  liteiature,  many  relatively 
common  problems  with  free  surfaces  and/or  cavities  have 
not  been  solved  in  closed  form  without  introducing  a 
number  of  simplifying  assumptions  which  are  not  in 
accord  with  real  situations.  Available  analytic  methods 
generally  require  that  the  fluid  be  assumed  weightless  (i.e. 
the  acceleration  of  gravity  is  zero).  Furthermore,  since 
such  methods  are  based  on  complex  variables,  they  are 
restricted  to  plane  two-dimensional  flows.  Consequently, 
in  order  to  solve  problems  with  free  surfaces  under  the 
influence  of  gravity,  or  three-dimensional  problems  even 
if  axially  symmetric,  researchers  have  been  forced  to 
obtain  solutions  based  on  numerical  approximations 
rather  than  solving  the  problems  in  closed  form. 

The  application  of  finite  differences  constitutes  one 
of  the  most  powerful  and  universally  applicaole  methods 
for  obtaining  such  approximate  solutions.  The  use  of 
finite  differences  for  solving  free  streamline  problems  in 
the  physical  plane  is  extremely  difficult  since  the  position 
of  the  free  streamlines  is  unknown  a  priori  The  solution 
can  be  obtained  only  through  a  process  of  repeatedly 
adjusting  the  assumed  position  of  the  free  streamlines, 
through  considerable  insight  and  judgment,  until  all 
conditions  of  the  problem  are  satisfied.  Since  the  means 
for  determining  whether  all  conditions  are  satisfied  is 
often  quite  insensitive  to  the  position  of  the  free 
streamlines,  it  is  difficult  to  determine  the  reliability  of 
the  resulting  approximate  solution,  and  consequently  the 
literature  contains  a  number  of  examples  where  sub¬ 
sequent  analyses  have  demonstrated  that  considerable 
error  resulted  because  of  an  incorrect  position  of  the  free 
streamline. 

An  approach  for  solvin’  two-dimensional  fluid  flow 
problems  which  is  supe.  ior  in  many  regards  to  a  formula¬ 
tion  in  the  physical  plane,  particularly  if  free  surfaces  are 
present,  is  to  interchange  the  usual  role  of  variables  in  the 
problem.  Such  inverse  formulations  have  used  the  poten¬ 
tial  function,  <{> ,  and  the  stream  function,  v  -  as  the 
independent  variables,  and  as  dependent  variables  such 
quantities  as:  (1)  the  magnitude  of  the  cartesian  coor¬ 
dinates  x  and  y,  (2)  the  angle  of  the  direction  of  flow,  3. 
and  the  logarithm  of  the  magnitude  of  the  velocity,  log 
|V|  ,  or  (3)  the  magnitudes  of  the  horizontal  and  vertical 


components  of  the  velocity,  u  and  v.  A  major  advantage 
to  such  an  inverse  formulation  is  that  fr: .  surfaces,  being 
streamlines  along  which  is  constant,  become  straight 
boundaries  in  the  ifnji  plane,  and  many  problems  are 
consequently  confined  within  rectangular  regions.  Also 
the  results  from  a  solution  are  in  an  ideal  form  for 
plotting  a  flownet  and  are  well  adapted  for  computing 
other  quantities  of  interest  concerning  the  flow. 

This  type  of  inverse  formulation,  accompanied  by  a 
subsequent  finite  difference  solution,  has  been  used  to 
study  a  variety  of  two-dimensional  free  streamline  prob- 
i.  ms  (Thom  and  Apelt,  1961;  Cassidy,  1965;  Markland, 
1965;  Jeppson,  1966  and  1969a).  The  same  approach  has 
been  used  to  solve  problems  dealing  with  plane  saturated 
porous  media  flow  with  phreatic  or  free  surfaces  (Jepp¬ 
son,  1968a,  b,  and  c),  and  unsaturated  moisture  move¬ 
ment  in  soils  (Jeppson  and  Nelson,  1970,  and  Jeppson  et 
al.,  1972).  The  same  approach  of  using  a  lormulation 
which  interchanges  the  usual  role  of  the  variables  with  an 
accompanying  numerical  solution  has  been  used  to  solve 
three-dimensional  problems  with  axial  symmetry.  In  these 
problems  the  magnitude  of  the  radial  and  axial  coor¬ 
dinates  r  and  z  are  made  dependent  in  the  plane  of  the 
potential  function  and  Stokes’  stream  function  (or 
logarithm  thereof).  (See  Jeppson,  1966;  Mackenroth  and 
Fisher,  1968;  Jeppson,  1968d,  1969b  and  1970.) 

The  work  reported  herein  extends  the  inverse 
formulation  technique  which  has  been  used  in  solving 
plane  and  axisymmetric  potential  fluid  flow  problems  to 
general  three-dimensional  potential  fluid  flow  problems 
and  demonstrates  the  applicability  of  the  methods  by 
obtaining  a  numerical  solution  to  the  three-dimensional 
flow  in  an  open  channel  past  a  strut.  While  this  problem  is 
a  very  simple  three-dimenM.mal  problem  for  which  a 
two-dimensional  analysis  (or  for  some  features  a  one¬ 
dimensional  analysis)  may  be  adequate,  it  does  include  the 
common  boundary  conditions  found  in  most  problems. 
Furthermore,  because  of  the  simplicity  of  the  problem, 
the  adequacy  or  inadequacies  of  the  numerical  solution 
can  more  readily  be  ascertained  and  where  necessary, 
modifications  made.  Conseouently  the  results  from  the 
problem  solution  have  the  primary  purpose  of  illustrating 
the  method  of  inversely  formulating  and  solving  a 
three-dimensional  free  surface  flow  problem.  With  a  better 
understanding  of  the  performance  of  various  numerical 
schemes  in  solving  inversely  formulated  three-dimensional 
problems,  the  next  step  would  be  to  apply  the  methods  to 
more  complex  three  dimensional  flows. 
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INVERSE  FORMULATION 


Selection  of  variables 


The  first  step  in  developing  an  inverse  formulation 
to  three-dimensional  potential  flows  is  the  selection  of 
three  appropriate  dependent  and  three  appropriate 
independent  variables.  Since  the  best  inverse  approaches 
in  the  literature  to  plane  and  axisymmetric  flow  problems 
have  considered  the  magnitudes  of  the  coordinates  x  and 
y  or  r  and  z  as  dependent  variables,  the  magnitudes  of  the 
cartesian  coordinates  x,  y,  and  z  should  constitute 
appropriate  dependent  variables  in  an  inverse  formulation 
to  a  three-dimensional  problem.  This  same  literature 
suggests  that  the  potential  function  as  well  as  some 
functions  to  define  the  flow  paths  would  constitute 
appropriate  independent  variables,  or  define  the  space 
within  which  the  problem  is  defined.  The  functions 
selected  for  defining  the  flow  paths  consist  of  two  stream 
surfaces  which  are  tangent  to  the  velocity  vector.  Yih 
(1957)  has  given  equations  for  defining  two  such  stream 
functions  which  will  be  denoted  by  *!>  and  in  this 
report.  Nelson  (1963)  has  given  equivalent  definitions  for 
use  in  three-dimensional  porous  media  flow  applications. 


The  basic  equations  in  these  definitions  are: 


=  *yC  = 

•  •  •  -  (I) 

II 

*s 

i* 

II 

.  -  .  -  (2) 

=  -  *7. 

.  .  .  .  (3) 

in  which  p.  v,  and  w  are  the  components  of  the  velocity 
vector  V  in  the  x.  and^  z.  coordinate  directions 
respectively,  i.e.  V  =  ui  +  vj  +  wT.  and  the  subscripts 
denote  partial  derivatives  in  the  usual  manner,  i.e.  '*r  = 
3H'  /  3z.  etc.  It  can  easily  be  shown  that  F.qs.  1.  2.  and  3 


reduce  to  the  well  known  equations 


V 


and  V  = 


$y  for  the  special  case  of  plane  potential  flows.  In 
vector  notation  Eqs.  1 . 2.  and  3  become 


V  =  (grad  *)  *  (grad  )  =  grad  S’ 


(4) 
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Transformation  from  physical 
space  to  Qyy  *  space 


To  obtain  the  basic  inverse  equations  note  that  since 
the  potential  function  and  the  two  stream  functions  are 
functions  of  x,  v,  and  z,  i.e.  $  =  F(x,y,z),  ¥  =  G(x,y,z) 
and  V*  =  H(x,y,z),  it  follows  that  x,  y,  and  z  must  also  be 
functions  of  $ ,  4\  and  'V*,  i.e.  x  =  f(  $,  4',  f*),  y  = 
g(  ❖.I’,  4*),  and  z  =  h(4>,4',  4*).  Using  the  chain  rule  to 
differentiate  x  =  f(  $.  4’,4'  *)  with  respect  to  x.  y,  and  z 
respectively  gives 


1  =  x*Fx  +  x$Gr  +  Xr*H 


X 


0  =  x  F  +  x,G  +  xt*H 

y  ^  y  i  y 


0  =  +  xt*H 


+  y*,  +  v 


(5) 


Solving  these  three  equations  for  the  unknowns  .  x^, 
and  x^.j.  gives 


-  -L  _  _  ±  d(F,H) 

J  d(y,  z)  ’  X$  J  d(y,  z)  ’ 


.  1  d(F.G) 

an  x-*  -  j  5(y7 z) 


(6) 


in  which  J  is  the  Jacobian  given  by  the  determinant 


J  = 


G  G  G 


H  H  H 


and  tin'  derivatives  of  the  quantities  enclosed  in  paren¬ 
theses  denote  minor  determinants  in  the  usual  wav.  i.e. 


'(G.  HI 


•(y.  r.) 


=  G  H  -G  It 

y  r.  7.  y 


Differentiating  v  =  g(  ;  .? .  »  *)  with  respect  to  \.  y. 
and  i.  respectively  and  solving  the  three  equations  gives. 


HaBUftt 


M 
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J H  w 11  ■  1  ■  ■' 


_l.  ajG.Hj  _  i  d(F,H) 
y4>  J  d(x.  z)  *  yi  "  J  d(x.  z)  ’ 


These  three  equations  are  the  basic  inverse  equations 
which  define  the  dependent  variables  x,  y,  and  z  in  the 
$4141*  space,  just  as  Eqs.  1,  2,  and  3  are  the  basic 


l  9(F,G)  equations  for  <t>,  ¥,  and  ¥*  in  the  physical  space. 

'  j  d{x,  z)  . (7)  Consequently,  when  associated  with  appropriate  bound- 


Likewise  differentiation  of  z  =  h( 4>, ¥,¥*)  leads  to. 


I  g(G,H)  1  &(F,H) 

J  5(x,  y)  *  *$  "J  d(x,  y)  * 


I  d(F,G) 
J  d(x,  y) 


Following  the  same  procedure  as  that  above  but 
solving  for  4  «,  4  y  ¥*  gives 


*  .  1  |jy.  *L  «  =  .  1  for,*)  * 

x  j  b($,$*)  y  j  5($7$*)  *  j  §($.$*) 

1  -  .  .1  JUXiiL '  q,  ,  4  -  .  1 

“x  j  0(575*)  ”y  j  57®.$*)’  -r  j  5(5,$*) 

a*  =  1  a*=  a*_  1  jfo.y) 

*  i  5(57*7  *  t  j  5(57$j '  ^  "  j  5(57$)  •  (9) 


in  which  j  is  the  inverse  Jacobian  determinant 
f®  f$  f$* 

j  *  8|»  =  J 

h<t>  h$  h^» 

By  substituting  from  Eqs.  6, 7, 8,  and  9  into  Eqs.  1 , 
2,  and  3,  the  following  three  inverse  equations  are 
obtained: 


xj>  =  y$  -  y^» 


y$  =  ■  *$*$» 


*♦  *  *$T$* "  *$»y| 


ary  conditions  for  a  particular  problem,  the  simultaneous 
solution  of  Eqs.  10,  11,  and  12  would  constitute  the 
solution  to  that  particular  problem.  Before  discussing 
methods  for  solving  these  equations  some  properties  of 
the  stream  surfaces  ¥  and  ¥*  will  be  presented. 

Properties  of  stream  surfaces 

The  definitions  for  stream  functions  ¥  and  HI*  as 
given  by  Eqs.  1,  2,  and  3  (or  Eq.  4)  satisfy  the 
incompressible,  steady  state  continuity  equation  V  *\T  = 
0.  This  can  be  verified  fromjhe  vector  identity  7*(A  x 
B}=  iTl  (7  x  A}-A.(V  x  B).  Thus  from  Eq.4, 


v.v  = 


;  V4-I  -  0$  .(V 


but  the  cur!  of  the  gradient  of  any  scalar  function  is  zero 
and  therefore  VxV1?  =  0  and  V  x  V¥*  =  0,  with  the 
result  that 


yv  =  v*(FJ  *v$*)  =  0 


The  stream  surfaces  defined  by  holding  both  HI  and 
HI  *  constant  are  orthogonal  to  the  equipotential  surfaces 
defined  bv  holding  $  constant.  Orthogonality  exists 
provided  the  dot  products  of  the  gradients  are  identically 
equal  to  zero.  Using  Eqs.  1,  2,  and  3  it  can  readily  be 
shown  that  V4»V¥  =  0  and  V4  *  V¥*  =  0  and 

therefore  the  equipotential  surfaces  are  everywhere  at 
right  angles  to  the  surfaces  defined  by  holding  the  two 
stream  functions  constant. 

In  general,  the  definitions  for  H/  and  HI*  do  not 
require  that  the  surfaces  defined  by  holding  ¥  and  ¥* 
constant  are  orthogonal  to  each  other.  However,  in  the 
previously  given  inverse  equations  it  is  necessary  that  of 
the  many  ¥  and  ¥*  equal  constant  surfaces  which  exist, 
only  those  are  selected  which  constitute  an  orthogonal 
pair  so  that  the  inverse  coordinates  4,  ¥,  and  ¥*  are 
independent.  The  use  of  the  inverse  formulation  assumes 
that  using  4,  ¥ ,  and  4*  as  orthogonal  coordinates  insures 
that  appropriate  orthogonal  ¥  and  ¥*  stream  surfaces  are 
selected.  Perhaps  a  more  fundamental  approach  would 
impose  the  condition  V  ¥*V¥*  =  0  directly.  Methods 
for  imposing  this  condition  directly  are  not  apparent, 
however. 
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METHODS  FOR  SOLVING  INVERSE  EQUATIONS 


Alternatives  available 

Considerable  guidance  in  the  selection  of  the  in¬ 
dependent  and  dependent  variables  for  the  three- 
dimensional  problem  was  provided  by  past  inverse  solu¬ 
tions  to  plane  flow  problems.  Less  guidance  is  available 
from  these  past  solutions  regarding  the  best  methods  for 
solving  the  three-dimensional  inverse  equations.  Since  the 
basic  equations  (Eqs.  10,  11,  and  12)  arc  nonlinear,  and 
each  equation  contains  all  three  dependent  variables  x,  y, 
and  z,  it  is  clear  that  numerical  methods  offer  the  best 
presently  available  approach  to  a  solution.  In  solving  the 


comparable  equations, 

**  =  y*  . (is) 

x*  ■  'yf  . (16) 


from  plane  potential  flows,  the  equations  are  first 
combined  by  differentiation  to  obtain  equations  involving 
only  one  dependent  variable.  These  equations  for  plane 
flows  are  the  inverse  Laplacian  equations  =  0  and 

V.  ?  y  =  0.  Because  of  the  products  present  in  the  terms 
onthe  right  side  of  the  equal  sign  in  Eqs.  10, 1 1 ,  and  12, 
reasonably  simple  equations  with  only  one  dependent 
variable  cannot  be  obtained  by  differentiation  and  com¬ 
bination  as  can  be  done  for  the  equivalent  plane  flow 
equations.  Consequently  an  alternative  approach  for 
solving  the  inverse  three-dimensional  equations  must  be 
sought. 

An  alternate  which  may  appear  feasible  at  first 
would  utilize  finite  difference  methods  to  obtain  a 
simultaneous  solution  to  the  boundary  value  problems 
associated  with  the  three  first  order  partial  differential 
Equations  10,  11,  and  12.  An  examination  of  the  finite 
difference  equations  obtained  from  these  three  equations 
by  approximating  the  derivatives  with  second  or  higher 
order  differences  indicates  that  point  by  point  iterative 
methods,  such  as  Gauss-Seidel  or  SOR  method  would  not 
be  convergent.  Such  iterative  methods  would  diverge 
simply  because  the  coefficient  associated  with  the  value  of 
the  variable  at  the  grid  point  in  question  is  less  in 
magnitude  than  the  sum  of  the  coefficients  of  the  other 
terms.  In  a  linear  system  the  equivalent  would  be  a 
nondiagonally  dominant  coefficient  matrix.  But  since 
diagonal  dominance  is  a  necessary  condition  for  point  by 
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point  methods  to  converge,  it  may  be  concluded  that  only 
if  first  order  forward  or  backward  differences  are  used  to 
replace  the  derivative  in  Eqs.  10,  11,  and  12  would  it  be 
possible  to  solve  the  boundary  value  problems  associated 
with  the  first  order  equations  simultaneously.  Because  of 
the  low  order  approximation  of  first  order  differences  this 
possibility  for  solution  was  not  considered  initially.  (Using 
a  weighting  all  possible  first  differences  which  depend 
upon  the  distance  from  the  boundary,  a  workable  method 
results.  This  approach  is  under  investigation  in  the  same 
project.)  Rather  three  alternatives  were  investigated. 
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The  first  alternative  is  to  use  block  iterative  meth¬ 
ods  to  solve  the  finite  difference  equations  obtained  from 
third  order  difference  approximations  of  the  derivatives  in 
the  first  order  partial  differential  equations.  The  merits  of 
this  approach  were  actually  investigated  by  implementing 
its  use  in  computer  programs  which  solved  the  finite 
difference  equations  across  an  entire  line  of  grid  points, 
and  across  two  adjacent  lines  simultaneously  for  the 
two-dimensional  problems  of  corner  flow.  The  conclusion 
was  that  these  block  (i.e.  line)  iterative  methods  were  also 
nonconvergent.  Later  study  has,  however,  shown  that 
what  was  considered  nonconvergence  may  have  actually 
been  due  to  the  poor  approximation  of  the  finite 
difference  solution  to  the  actual  corner  flow.  Regardless 
of  the  incorrectness  of  the  above  conclusion,  the  use  of 
olock  iterative  methods  for  solving  the  simultaneous 
boundary  value  problems  was  not  pursued  further.  Rather 
the  method  of  approach  which  is  described  in  this  report 
was  developed  and  implemented  in  a  computer  program 
for  solving  three-dimensional  flow  around  a  strut. 
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The  third  alternative  which  has  been  studied  for  4 

solving  the  simultaneous  boundary  value  problems  will  be  » 

described  more  fully  in  a  subsequent  report  containing  the  | 

results  of  a  Ph.D.  thesis  by  Allen  Davis.  Basically  this  later  '• 

alternative  also  uses  Eqs.  10,  11.  and  12  in  their  present  J 

forms,  and  obtains  a  simultaneous  solution  for  x.  y.  and  z  J 

from  the  difference  equations  obtained  by  third  order  | 

approxin.  itions  at  all  grid  points  on  an  entire  plane  within  | 

the  flow  space.  The  resulting  finite  difference  equations  | 

become  linear  under  the  assumption  that  values  on  | 

adjacent  planes  arc  known.  Consequently  the  solution  on  | 

each  plane  can  be  obtained  efficiently  by  utilizing  1 

techniques  for  grouping  tire  m  nzero  elements  of  the  | 

coefficient  matrix  into  hands,  and  implement  iifg  | 


^i*,*****^^^  .  ■■  ...  ■  , . ,:4,..,..-^.. 


the  matrix.  By  repeatedly  obtaining  such  solutions,  plane 
by  plane  and  subsequently  repeating  the  entire  process 
until  the  absolute  sum  of  changes  in  all  three  variables  at 
all  grid  points  became  less  than  some  error  parameter,  the 
final  solution  should  result.  In  essence  this  alternative  is 
an  extension  of  block  iterative  methods  to  a  space 
boundary  value  problem  in  which  the  block  becomes  an 
entire  plane  and  direct  methods  are  used  to  solve  the 
finite  difference  equations  in  that  plane.  For  some  yet 
unknown  reason  this  procedure  neither  converges  to,  nor 
diverges  from,  the  final  solution.  A  more  detailed  descrip¬ 
tion  of  the  implementation  of  this  method  and  its 
inadequacy  will  be  given  in  a  subsequent  report. 

The  method  of  solution  described  in  this  report 
does  not  use  Eqs.  10,  11,  and  12  directly  in  their  present 
form.  Rather  these  three  equations  are  combined  by 
differentiation,  under  the  assumption  that  certain  of  the 
derivatives  are  known,  in  such  a  way  that  quasi-separate 
equations  are  developed  for  each  variable  x,  y,  and  z  in 
different  planes  within  the  space.  The  magnitude 

of  the  assumed  known  quantities  in  these  separate 
equations  can  only  be  determined  approximately  until  the 
final  solution  is  obtained.  Consequently  these  assumed 
known  quantities  are  repeatedly  adjusted  in  a  cycle  of 
solutions  until  their  correct  values  are  obtained. 


a  d  e  f  g 
/  /  /  /  / 


cl\> 

.  .  .(20) 

b 

h  g  i  e 

/ 

/  /  /  / 

V* 

=  3V,Z+  ■  x»(iz+*  •  • 

.  .  .(21) 

c 

if  h  d 

/ 

/  /  /  f 

.  .  .  .(22) 

clz« 

=  x+v  -  •  • 

in  which  c,  =  D/M,  ,  and  the  single  letter  over  the 
individual  derivatives  will  be  used  subsequently  whenever 
that  derivative  is  assumed  to  be  known. 

Development  of  quasi-separate 
equations  for  x,  y,  and  z 

To  demonstrate  how  separate  quasi-separate  equa¬ 
tions  for  x,  y,  and  z,  which  apply  on  an  individual  plane 
within  the  *  space,  might  be  developed,  Eqs.  20  and 
21  are  written  below  assuming  that  derivatives  with 
respect  to  rj>*  are  known  and  that  the  variable  z  is  known. 


Nondimensionalizing  independent  variables 

Before  demonstrating  how  such  quasi-separate  equa¬ 
tions  can  be  obtained,  Eqs.  10,  11,  and  12  will  be 
transformed  so  the  new  independent  variables  <p,  and 
V*  are  dimensionless  as  given  by  the  following  three 
equations 


cix*=eVfg . (20a) 

cl*$  *  '  ex+ . (21a) 


Upon  differentiating  Eq.  20a  with  respect  to  $  and 
Eq.  21a  vwth  respect  to  $  and  elimi'ialmg  y  =  y 
gives  W  ^ 


= 


+  = 


N,D4> 


M1  $ 


.(17) 


n/Q 


.(18) 


c.x**+  7-  +  -  --  <C,'V  fg) 

. (23) 


'  Cj  [  *44  ,0'  r+ 

+  (fg)d  =  0 


sS 


-i 


•f  = 


gig 

slQ 


Likewise  differentiating  Eq.  20a  with  respect  to  ij> 
and  Eq.  21a  with  respect  to  <J>  and  combining  the 
0  ■’)  resulting  equations  gives  the  following  equation  for  y  in 
planes, 


i 

a 


3 

"S 


in  which  D  is  the  undisturbed  depth  of  flow  in  the 
channel,  Q  is  the  total  volumetric  flow  rate  in  the  channel, 
N,  is  the  number  of  grid  increments  in  the  <J>*  direction 
and  M,  is  the  number  of  grid  increments  in  the 
direction.  Transforming  Eqs.  10, 11,  and  12  by  means  of 
Eqs.  17, 18,  and  19  leads  to 


ciy*a  +  f;  V+  *  (fg)+  ]  +  T(gh‘ ciV 


-  <gl»*  «  0 


(24) 
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If  y  had  been  considered  known  along  with  deriva¬ 
tives  with  respect  to  ip*  in  Eqs.  20  and  22  the  same 
procedure  as  that  given  above  would  have  resulted  in 
equations  for  x  and  z  in  the  <f>ip  planes.  Equations  for  y 
and  z  would  result  if  x  were  considered  known.  In  fact  for 
each  pair  of  equations  which  can  be  formed  from  Eqs.  20, 
21,  and  22,  two  quasi-separate  equations  would  result 
under  each  assumption  of  known  variables.  Table  1  lists 
the  18  equations  that  can  be  obtained  in  this  manner,  and 
indicates  in  which  plane  each  of  these  equations  applies  as 
well  as  from  which  two  of  the  three  basic  equations  it  was 
obtained. 

The  motivation  for  combining  Eqs.  20,  21,  and  22 
by  differentiation  into  the  second  order  partial  differen¬ 
tial  equations  in  Table  1  is  to  obtain  second  derivatives  in 
the  equation,  for  which  second  order  central  difference 
approximations  lead  to  diagonally  dominant  coefficient 
matrices.  The  equations  in  Table  1  also  have  some 
resemblance  to  Laplacian  type  equations  for  which  the 
common  finite  difference  methods  have  been  developed. 
Perhaps  the  greatest  motivation  for  developing  the  equa¬ 
tions  in  Table  1 ,  however,  was  to  have  separate  equations 
from  which  to  solve  each  of  the  dependent  variables  x,  y, 
andz. 

Criteria  for  selecting 
best  suited  equations 

Generally  in  solving  any  particular  problem,  only 
one  equation  for  each  of  the  unknowns  x.  y,  and  z  would 
be  used.  Should  considerable  differences  exist  in  the  flow 
patterns  in  different  regions  of  a  particular  problem  it 
may  be  desirable  to  use  different  equations  in  different 
regions.  The  success  and  efficiency  of  obtaining  a  solution 
by  use  of  the  equations  in  Table  1  depends  upon  the 
selection  of  the  equation  which  will  be  used  to  solve  for 
each  of  the  unknowns.  While  there  are  additional  criteria 
which  might  help  in  making  this  selection  it  appears  that 
the  following  three  items  are  important:  (1)  The  assump¬ 
tions  of  known  derivatives  should  be  made  as  valid  as 
possible;  that  is  the  values,  denoted  by  single  letters  in  the 
equation  that  is  used,  should  be  maintained  as  constant  as 
possible  during  the  solution  process  which  would  start 
with  an  initialization  and  proceed  until  all  conditions  of 
the  problem  are  satisfied.  (2)  That  the  coefficients  of  the 
two  second  derivative  terms  in  the  equation  be  as  nearly 
equal  as  possible.  For  the  last  three  equations  in  Table  1. 
obviously  at  least  one  of  the  single  lettered  values  must  be 
negative  so  that  the  PDE  is  elliptic.  (3)  That  the 
magnitudes  of  the  terms  involving  first  derivatives  be 
maintained  as  small  as  possible. 

Several  reasons  exist  for  citing  these  criteria.  First  if 
the  single  lettered  values,  which  are  assumed  to  be  known 
during  the  process  of  obtaining  a  solution  on  any  plane, 
have  their  values  altered  greatly  between  successive 
solutions  in  that  plane,  they  will  obviously  affect  the 
results  from  these  consecutive  solutions.  These  solution 
results,  in  turn,  could  affect  the  magnitudes  of  the 


coefficients.  If,  however,  the  lettered  coefficients  are 
nearly  correct  at  initialization,  or  if  their  magnitudes  have 
a  minor  influence  on  the  solution,  the  final  solution  will 
be  obtained  in  fewer  total  arithmetic  operations. 

The  basis  for  the  second  criteria  is  to  make  the 
solution  in  each  of  the  planes  of  the  space  equally 
dependent  upon  all  four  of  its  boundary  values  (or 
boundary  conditions),  and  not  more  dependent  on  two 
opposite  boundaries  than  on  the  other  two  opposite 
boundaries.  This  latter  condition  would  occur  if  the 
coefficient  of  one  of  the  second  derivatives  was  very  small 
in  comparison  to  the  other.  The  second  criteria  also  helps 
insure  that  the  equation  has  some  resemblance  to  Laplaces 
equation  for  which  many  numerical,  as  well  as  other, 
solutions  have  been  obtained.  Should  this  criteria  be 
strongly  violated,  a  solution  in  each  individual  plane  may 
be  obtained  with  fewer  numerical  calculations  by  simul¬ 
taneously  solving  the  system  of  finite  difference  equations 
along  the  grid  lines  in  the  direction  of  the  independent 
variable  whose  derivative  has  the  larger  coefficient.  Since 
the  problem  is  of  the  elliptic  type,  this  would  mean  that  a 
high  dependency  must  exist  between  the  values  on  this 
plane  and  those  on  adjace  t  planes.  Consequently,  any 
reduction  in  arithmetic  calculations  in  obtaining  individ¬ 
ual  solutions  would  be  more  than  offset  by  more  cycles  of 
such  solutions.  Furthermore,  the  solution  process  may  be 
less  likely  to  be  convergent.  Consequently  satisfying  the 
second  criteria  simultaneously  helps  assume  that  the  first 
criteria  is  satisfied. 

An  illustration  of  how  these  criteria  aid  in  the 
selection  of  the  equation  which  will  be  used  to  solve  each 
of  the  dependent  variables  x,  y,  and  z  is  given  later  in  the 
discussion  of  the  problem  of  flow  around  a  strut  in  a 
channel. 

To  obtain  a  better  understanding  of  how  rapidly,  or 
whether,  iterative  finite  difference  methods  would  con¬ 
verge  for  the  equations  in  Table  I,  individual  computer 
programs  were  written  to  solve  each  of  the  equations  in 
Table  1.  For  each  such  problem  Diridilet  boundary 
conditions  were  specified,  and  algorithms  implementing 
both  the  successive  over-relaxation  (SCR)  method  and  the 
line  successive  over-relaxation  (LSOR)  method  (see 
Forsythe  and  Wasow,  1960.  or  Vaiga,  1962)  were  tested. 
In  part  the  criteria  given  lor  selecting  the  best  equations 
for  solving  a  particular  problem  were  arrived  at  from 
noting  and  comparing  the  performance  of  these  separate 
programs  in  obtaining  a  solution.  The  performance  of 
those  programs  implementing  solutions  to  Eqs.  35 
through  40  in  rv*  planes  was  generally  considerably 
poorer  than  those  implementing  solutions  in  either  +i>  or 
S  v*  planes.  If  a  poor  initialization  of  all  unknowns  was 
used  when  solving  the  equations  in  d't^*  planes,  solutions 
did  not  result,  but  rather  rapid  divergence  occurred.  When 
such  lack  of  convergence  occurred  it  appeared  *o  be 
associated  with  initializations  which  at  some  grid  points 
caused  the  coefficients  of  the  second  derivative  terms  to 
have  opposite  signs,  or  which  caused  the  magnitudes  of 
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Table  1.  Quasi -separate  equations  obtained  from  Eqs.  20, 21 ,  and  22  by  assuming  some  variables  were  known. 


Derived 

from 

Eqs. 

Plane 

of 

Equation 

23 

20  It  21 

44> 

24 

20  It  21 

44> 

25 

20  It  22 

44* 

26 

20  It  22 

27 

21  It  22 

*4- 

28 

21  It  22 

0+ 

29 

20  fc  21 

*>* 

30 

20  It  21 

♦+* 

31 

20  It  22 

44* 

32 

20  It  22 

H* 

33 

21  It  22 

♦4>* 

34 

21  It  22 

♦+* 

35 

20  It  21 

44* 

36 

20  It  21 

44* 

37 

20  It  22 

44* 

38 

20  It  22 

44* 

39 

21  It  22 

44* 

40 

21  It  22 

+4* 

Quasi -Separate  Partial  Diff.  Eq. 

cix<M  +  7[  [“V  V+  ‘  <ghM  -  T-  <ci *♦ + +  <fg)*  =  0 
ciyW + 7[  teyw+  V* '  (fg)+]  +  *T  (gb  *  ciV  - <gh>*  *  0 
cix**+i;  <%*+W  (dh,+1+  i (de  -ciV  -(de)*=  0 
v«+i;Pv+,A-,deVT  (V#+dh,  +  (dh,^  o 

ciy** thyH+ V*  * (i£ V  *  TT  (ciy*  +  ie)  +  (le)*  =  0 

C1*W+^  f“V  V*‘  <ie)*>  +  IT  (lf  -ciV  *  (in*  =  0 
cixoe+  f[  -**+•♦* +  VV  ‘  <io)+*1+ if (de '  ciV  ‘ (de>*  *  0 
ciW^  t*W*+  vv,dV^'  (ciV ie)  +  (ie)*  =  0 

.  d 

etx**+^  l*1*****  V**  ‘  (iV  +  d  (ciV  fg,+  W*  *  o 
ci**+^  <dV+*+  w  *  (,g,+*1+'d  w*ciV  +  (if)e  =  0 
ciy^4+—  +  W  *  <dh)+*]+T  (gh~ciV‘(,h,o  =  0 

V**f  P*w  +  W  ■  (ghVJ'  *?  (ciV  dh,+  (dlV  0 

a 

tha^  -  a*w  +  (dh+„  +  hf+>x+-  <di+„  +  hd+)a+,  -  Cj(db+*  +  ha+)  *  0 

“‘V^V***  (%+  VN ■  (iV+  W  +  C1(V‘%)  * 0 

«brH-  rr^*  *■  <*%+  «fr+*>r+  -  <*t+*+ »»*+>r+*+  0 
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these  coefficients  to  take  on  small  values.  From  this 
experience,  the  additional  criteria  may  be  added  to  the 
above  three  to  avoid  in  general,  one  of  the  equations 
which  apply  on  planes. 

The  performance  of  the  two  algorithms  (SOR  and 
LSOR  methods)  which  were  implemented  in  the  programs 
mentioned  above  indicated  that  generally  less  computer 
execution  time  was  required  by  the  LSOR  method  than 
the  SOR  method.  No  experimentation  was  done  to 
examine  the  effect  of  the  over-relaxation  factor.  All 
solutions  under  both  methods  used  an  over-relaxation 
factor  equal  to  1.4.  Since  the  comparison  was  close, 
showing  that  the  LSOR  method  required  in  the  neighbor¬ 
hood  of  20  to  40  percent  less  execution  time  than  the 
SOR  method,  changes  in  over-relaxation  factors  might 
favor  the  SOR  method.  Furthermore,  the  outcome  of 
such  a  comparison  is  computer  system  dependent,  as  well 
as  being  influenced  by  the  particular  statements  in  the 
source  language  written  for  each  method.  In  the  LSOR 
method  more  computations  are  involved  per  iteration  but 
fewer  iterations  are  required  for  a  solution  than  with  the 
SOR  method.  Since  the  additional  computations  per 
iteration  are  primarily  with  nonsubscripted  variables  or 
single  arrayed  variables,  the  LSOR  method  requires  fewer 
operations  with  triple  subscripted  arrays.  These  compari¬ 
sons  were  made  on  the  UNIVAC  1108  system,  under 
EXEC  8,  at  the  University  of  Utah. 

Finite  differences 

The  finite  difference  operators  have  been  obtained 
by  replacing  the  derivatives  in  the  equations  in  Table  1 
with  second  order  central  difference  approximations.  The 
finite  difference  space  network  has  been  selected  such 
that  AiJ;*  =  Aiji  =  A<J>  =  1.  The  grid  spacing  increments 
A  -,]>*,  A’i',  and  A$  can  each  be  equated  to  unity  because 
the  number  of  grid  increments.  M,  and  N,  in  the  ^  and 
•£*  directions,  respectively,  are  included  in  Eqs.  17,  18, 
and  19  for  defining  the  dimensionless  coordinates  $  , 
and  4-*.  The  motivation  for  introducing  M,  and  N,  in 
Eqs.  17,  18,  and  19  was  to  allow  these  increments  to  be 
unity  and  thus  eliminate  a  number  of  multiplications 
which  would  result  from  nonunity  A’s  in  the  finite 
difference  operators. 

The  finite  difference  operators  for  interior  space 
grid  points  are  given  in  Table  2  for  the  first  12  equations 
in  Table  1.  To  make  for  easy  reference  the  equation 
numbers  given  for  each  finite  difference  operator  in  Table 
2  is  the  same  as  that  for  the  PDE  from  which  it  was 
derived  in  Table  1 .  The  forms  of  these  finite  difference 
operators,  as  given  in  Table  2,  conform  to  that  needed  to 
apply  the  LSOR  method  along  those  lines  defined  by  the 
incremented  subscripts  on  the  left  side  of  the  equal  sign 
(i.e.  in  the  direction  of  increasing  )  and  which  lie  in  the 
plane  on  which  the  particular  equation  applies.  The  triple 
subscripts  to  x,  y,  and  z  in  the  finite  difference  operators 
are  arranged  in  order  so  that  the  first  corresponds  to  4. 
the  second  to  ^  and  the  third  to  f*  as  defined  by 


i  =  1  +  4/A0  (41) 

j  =  l  +  +/A9  (42) 

k  =  1  +  +S/A4? . (43) 


The  oi’s  in  each  finite  difference  operator  are  unique  to 
that  operator  as  defined  in  the  right  portion  of  the  table. 
They  are  used  to  simplify  the  writing  of  the  operators  and 
consist  of  the  combination,  and/or  derivatives,  of  the 
assumed  temporarily  known  quantities  given  by  a  single 
letter  in  the  PDE’s  in  Table  1. 

The  operators  in  Table  2  can  be  rewritten  readily  to 
conform  to  that  needed  tc  apply  them  in  the  LSOR 
method  in  the  other  coordinate  direction  by  interchanging 
the  terms  across  the  equal  signs  or  for  the  SOR  method  by 
placing  only  the  term  with  an  ijk  subscript  on  the  left  of 
the  equal  sign. 

Numerical  operations  involved 
in  obtaining  a  solution 

As  pointed  out  earlier,  the  solution  or  solutions  on 
any  given  plane  within  the  space  of  the  problem  must  be 
obtained  repeatedly;  each  subsequent  solution  hopefully 
will  have  more  nearly  correct  coefficients  which  are 
assumed  known,  but  which  actually  are  dependent  upon 
knowing  the  correct  solutions  to  the  other  variables. 
Consequently,  a  single  group  of  solutions  on  individual 
planes  for  x,  y.  and  z  will  not  be  sufficient.  Rather  such 
groups  of  solutions  must  be  obtained  repeatedly  until  all 
coefficients  are  correct.  To  help  describe  the  procedure 
used  in  obtaining  the  final  complete  finite  difference 
solution  to  a  three-dimensional  problem,  the  following 
terminology  will  be  used. 

(a)  Tentative  solution -teieis  to  a  solution  based  on 

any  of  the  finite  difference  operators  in  Table  2  (or  any  of 
the  finite  difference  operators  for  a  boundary  condition  as 
given  later)  on  a  specified  plane  within  the  space. 

These  tentative  solutions  3re  based  on  the  best  values  for 
the  ex's  which  can  be  computed  at  that  stage  in  the 
solution  process,  and  consequently  they  are  only  tenta¬ 
tive.  but  the  results  of  these  solutions  are  needed  to 
obtain  better  estimates  of  the  os  in  the  operators  for  the 
other  two  dependent  variables. 

(b)  Iteration  number  refers  to  the  number  of  times 
the  LSOR-nethod  (or  whitever  other  method  is  being 
used)  adjusts  al<  the  values  of  x,  y.  or  z  on  the  particular 
plane  for  which  a  given  tentative  solution  is  being 
obtained.  A  sufficient  number  of  iterations  are  required 
for  each  tentative  solution  until  the  sum  across  all  grid 
points  of  that  plane  of  absolute  changes  in  value  of  that 
dependent  variable  is  less  than  the  prescribed  error 
parameter. 

(c)  Cycle  number  refers  to  the  number  of  times  all 
tentative  solutions  arc  obtained.  Thus  during  the  first 
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Table  2.  Finite  difference  operators  which  are  based  on  the  partial  differential  equations  in  Table  1. 


Eq. 

No. 

Finite  Difference  Operator 

Definition  of  a 

Coefficient*  in  Operator 
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Table  2.  Continued. 


Finite  Difference  Operator 
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cycle  all  tentative  solutions  for  x,  y,  and  z  will  be 
obtained  as  well  as  possibly  tentative  solutions  for  these 
variables  on  the  boundary  planes  which  are  not  of  the 
Dirichlet  type.  The  same  process  is  repeated  for  the 
second  cycle,  etc.  In  the  actual  computer  program  as  it 
has  been  written  for  the  problem  of  flow  around  a  strut, 
the  additional  capability  has  been  provided  to  repeat  all 
the  tentative  solutions  for  x,  y,  or  z  on  interior  planes 
more  than  one  time  before  proceeding  to  the  tentative 
solutions  of  the  next  variable. 

The  LSOR  method  has  been  used  for  the  reasons 
given  earlier  to  obtain  all  tentative  solutions  except  for 
certain  boundary  planes  as  will  be  described  later  for  the 
example  problem  which  is  given  herein.  While  a  descrip¬ 
tion  of  the  LSOR  method  can  be  found  in  a  number  of 
text  books  dealing  with  finite  difference  methods  for  PDE 
(see  for  example  ,ones,  1965),  a  brlrf  explanation  is  given 
here  for  the  sake  of  completeness  as  weU  as  to  point  out 
certain  unique  features  of  the  LSOR  method  as  applied  to 
the  equations  of  Table  2. 


In  implementing  the  algorithm  given  by  Eqs.  45  and  46  it 
is  not  necessary  to  set  aside  storage  for  a  new  array  f. 
Rather,  since  the  values  of  r  need  not  be  retained,  the 
values  of  f  may  be  stored  in  the  former  array  positions  for 


Upon  obtaining  the  solution  vector  X  which  repre¬ 
sents  the  unknown  values  across  an  entire  grid  line,  the 
individual  elements  are  immediately  adjusted  by  the 
over-relaxation  formula 


=  x.  +  w,(x.  -  xf., ) 


-  -  -(47) 


The  LSOR  method  can  be  understood  by  noting 
that  the  application  of  any  of  the  operators  in  Table  2 
across  all  interior  grid  points  of  a  line  leads  to  a  system  of 
linear  algebraic  equations  provided  that  the  cfs  and  values 
on  adjacent  lines  are  known.  In  matrix  notation  this 
system  has  the  form 


Ax  =  B 


(44) 


in  which  X  represents  the  unknown  vector.TTis  the  vector 
of  known  quantities  and  A  is  a  tridiagonal  matrix.  The 
fact  that  A  is  Iridiagcnal  is  an  important  feature  of  the 
method  from  a  computational  viewpoint,  since  such 
tridiagonal  systems  of  equations  can  be  solved  by  a  single 
pass  through  the  rows  with  a  Guassian  elimination  which 
leaves  a  matrix  with  only  two  elements  on  each  row;  the 
diagonal  element  and  the  next  element.  The  solution  to 
the  system  is  subsequently  obtained  by  back  substitution. 
Some  writers  have  referred  to  the  method  for  accomplish¬ 
ing  this  solution  as  the  Thomas  algorithm.  This  method 
defines  the  sequence  of  elements  of  A  immediately  to  the 
left  of,  on,  and  immediately  to  the  right  of  the  diagonal 
by  vectors  "if  77  and  Tjesgectively.  Then  additional 
elements  of  other  vectorsTand  g  are  defined  by 


t 

m 


- 712 - ,  R  =  f  (b  -  q  g  ,)/s 

r  -f  ,q  m  m  m  mm-!  m 
m  m-1  m 


2  <  m  s  n  .  .(45) 


in  which  n  is  the  number  of  rows  and  columns  in  A,  and 
the  b’s  are  the  elements  of  B?  The  solution  vector  X  is 
obtained  from 


in  which  the  Xj’s  (with  the  single  subscript)  are  the 
elements  of  IT,  and  x^k  (with  the  triple  subscript)  are  the 
values  of  the  dependent  variable  x,  y,  or  z  at  the  grid 
points  along  the  line  in  question.  The  superscript  p 
represents  the  number  of  the  iteration  and  w,  is  the 
over-relaxation  factor  with  a  value  between  zero  and 
unity.  Eq.  47  is  not  the  usual  form  of  the  over-relaxation 
equation  which  utilizes  an  over-relaxation  factor  w  =  w,+ 
1.  The  use  of  Eq.  47  in  place  of  the  more  usual  form  has 
the  advantage  that  the  computer  core  positions  for  the 
triple  subscripted  array  needs  to  be  located  once 
instead  of  twice  to  carry  out  the  arithmetic  on  the  right  of 
the  equal  sign. 

The  LSOR  method  proceeds  from  line  to  line  until 
the  value  of  the  variables  across  all  iines  within  the  plane 
have  been  adjusted.  Upon  completing  the  last  line  the 
entire  process  is  repeated  as  the  next  iteration,  etc.  In 
implementing  the  LSOR  method  for  the  equations  in 
Table  2,  the  as  have  been  computed  only  during  the  first 
iteration  and  stored  for  subsequent  iterations.  The  reason 
for  doing  this  becomes  obvious  upon  noting  that  some  of 
the  a’s  are  independent  of  the  solution  on  that  plane,  and 
those  that  are  have  minor  effect  on  the  resulting  solution. 
Consequently,  the  majority  of  the  arithmetic  involved  in 
solving  the  tridiagonal  system  repeatedly  is  with  single  real 
variables  or  single  subscripted  arrays.  Since  the  ds  will 
take  on  different  values  during  the  next  cycle,  particularly 
during  the  first  few  cycles  of  the  solution  process,  there  is 
no  need  to  iterate  ,u..'ii  the  tentative  solution  during  first 
cycles  satisfy  a  small  error  parameter.  By  permitting  a 
limited  number  of  iterations  to  occur  in  obtaining  any 
tentative  solution,  the  tentative  solutions  will  progessively 
satisfy  a  smaller  error  parameter,  until  eventually  during 
later  cycles  the  specified  error  parameter  will  be  satisfied 
with  fewer  than  the  maximum  specified  number  of 
iterations. 


FREE  SURFACE  FLOW  AROUND  A  STRUT 


In  the  initial  application  of  the  previously  described 
inverse  solution  to  three-dimensional  potential  flows,  a 
simple  problem  was  selected  to  test  the  workability  of  the 
methods.  The  first  such  problem  consisted  of  uniform 
flow  in  an  open  channel.  After  it  was  demonstrated  that 
the  methods  did  converge  to  the  solution,  providing  a 
reasonable  initialization  was  supplied  to  begin  the  process, 
the  program  used  for  a  solution  to  this  first  test  problem 
was  modified  to  solve  the  problem  of  flow  in  a  channel 
with  vertical  side  walls  past  a  vertical  strut  which  extends 
through  the  free  surface.  This  problem  is  described  in 
detail  herein.  It  demonstrates  possible  methods  for  in¬ 
corporating  boundary  conditions  into  the  solution  of 
three-dimensional  inversely  formulated  problems.  While 
the  problem  represents  what  one  might  refer  to  as  a 
“mildly  three-dimensional  problem,"  it  does  contain 
examples  of  the  commonly  encountered  boundary  condi¬ 
tions.  Besides  having  the  advantage  that  much  of  the  flow 
behavior  might  be  predicted  from  more  elementary 
analyses,  and  therefore  an  indication  of  the  adequacy  of 
the  methods  may  be  evaluated,  a  “mildly  three- 
dimensional  problem"  of  this  nature  provides  a  base  upon 
which  a  number  of  techniques  for  handling  different 
boundary  conditions  can  be  experimented  with  and  the 
best  of  the  alternatives  selected.  It  soon  became  apparent 
even  while  experimenting  with  the  first  problem  of 
uniform  channel  flow,  that  completely  satisfactory  meth¬ 
ods  for  handling  free  surfaces  or  cavity  surfaces  under  the 
influence  of  gravity  would  be  hard  to  come  by.  Hopefully, 
future  research  will  improve  upon  some  of  the  techniques 
described  herein. 

Formulation  and  boundary  conditions 

A  sketch  of  the  problem  of  channel  flow  past  a  strut 
in  the  physical  space  is  given  in  the  upper  portion  of  Fig. 
1  and  the  same  region  in  the  <j>iptp*  space  Is  given  in  the 
lower  portion  of  this  figure.  The  *  space  has  been 
selected  such  that  the  bottom  of  the  channel  defines  the 
t*i  =  0  stream  surface  and  the  top  free  surface  of  flow  is 
defined  by  the  Qip*  plane  obtained  by  holding  ^  =  M ,  = 
M-l  where  M,  is  the  numbe<  of  grid  increments  used  in 
the  finite  difference  solution  in  the  ip  direction.  (Remem¬ 
ber  Ai^  =  1.)  The  plane  defined  by  ip*  =  0 
corresponds  to  the  left  wall  of  the  channel  (when  facing 
downstream)  and  the  right  wall  by  iji*  =  N,  =  N-l.  The 
beginning  of  the  space  boundary  value  problem  through 
which  the  flow  enters  is  assumed  to  be  far  enough 


upstream  to  be  influenced  insignificantly  by  the  presence 
of  the  strut.  This  plane  is  given  by  4>  =  0,  and  the 
last  tyty*  plane  of  the  pipip*  boundary  value  problem  is 
defined  by  <j>  =  L ,  =  L-l . 

After  placing  the  problem  in  the  <J>M  *  space,  the 
next  step  in  the  formulation  consists  of  selecting  the 
equations,  from  those  in  Table  1,  to  be  used  to  solve  for 
each  of  the  three  dependent  variables  x,  y,  and  z.  This 
selection  should  be  based  on  consideration  of  the  criteria 
given  earlier.  For  the  problem  being  considered  here,  the 
major  component  of  velocity  throughout  all  except  small 
regions  near  the  front  and  rear  of  the  strut,  is  in  the 
x-direction  in  the  physical  space.  Furthermore  because  of 
the  placing  of  the  problem  in  the  space,  ihe 

channel  bottom  with  ip  =  0  is  normal  to  the  y-direction 
and  the  sides  of  the  channel  with  tp*  held  constant  are 
normal  to  the  z-direction.  Consequently,  greater  variation 
of  x  occurs  in  ipip  or  $i|/  *  planes,  than  in  planes. 
The  major  change  in  y  is  in  the  ^-direction.  Therefore,  a 
plane  defined  by  ip  as  one  coordinate  should  be  selected 
for  obtaining  the  solutions  for  y.  Likewise,  the  major 
variation  of  z  is  in  the  direction  of  if>*,  and  consequently 
z  would  be  fairly  constant  on  separate  <j>  >p  planes. 

Therefore,  the  first  criteria  stated  earlier,  namely 
that  the  assumption  of  knowns  be  as  valid  as  possible, 
would  suggest  that  x  could  be  solved  for  on  separate  4><l’ 
or  planes,  but  not  \p  \p  *  planes.  Clearly  the 

magnitude  of  x  ^  is  larger  than  or  in  general  and 
consequently  an  easily  generated  initialization  of  the 
problem  would  have  larger  errors  in  the  magnitudes  of 
x  <j>>  than  x^  orx^*. 

The  second  criteria,  namely  that  the  coefficients  of 
the  two  second  derivatives  have  nearly  equal  magnitudes, 
will  be  used  to  narrow  the  choice  down  further.  For  x  on 
pip  or  <W*  planes  the  available  equations  are  23, 25, 29, 
and  31.  In  comparing  the  magnitudes  of  coefficients  of 
second  derivatives  c,  may  be  compared  to  the  square  of 
the  single  letters  representing  the  derivative,  i.e.  with  c2 , 
f  *  g2,  and  d2.  Should  the  problem  be  specified  so  that 
the  magnitude  of  c,  is  close  to  unity,  as  will  be  the  case, 
then  either  Eq.  23  or  31  could  be  selected.  Equations  25 
and  29  are  eliminated  because  the  coefficients  f  =  fy  * 
and  g  =  Zjj,  are  much  smaller  whereas  e  =  z^*  ,  and  d  = 
yy  are  close  to  unity  in  magnitude.  The  final  selection 
between  Eq.  23  or  31  is  arbitrary.  In  solving  the  problem, 
Eq.  23  has  been  used. 
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The  equation  to  use  in  solving  for  y  is  between  24 
and  27  in  M  planes.  Since  the  magnitude  of  h  =  ju,*  is 
much  smaller  in  general  for  this  problem  than  c,,  Eq.  24 
will  be  used  to  obtain  the  tentative  solutions  for  y. 

For  obtaining  solutions  for  z  only  the  two  equations 
30  and  34  on  planes  will  be  considered  further.  Of 
these  two  equations,  34  is  eliminated  since  the  coefficient 
i  =  xq;  is  small  and  the  coefficient  d  =  y^  in  Eq.  30  is 
close  to  unity.  In  summary  the  tentative  solutions  for  x 
and  y  will  be  obtained  from  separate  #  planes.  Equation 
23  will  be  used  for  x  and  Eq.  24  for  y.  The  tentative 
solutions  for  z  will  be  by  means  of  Eq.  32  on  <p\p*  planes. 

Using  the  above  description  ~f  the  <$  \p*  space  and 
the  selected  equations,  the  boi.idarv  conditions  given 
below  can  be  developed.  Some  of  these  conditions  are  also 
shown  in  the  lower  portion  of  Fig.  1  to  help  identify  that 
boundary  and  its  condition  in  the  $  iJmJ^  space.  A 
description  of  how  these  boundary  condition  equations 
are  obtained  follows  the  listing  of  the  equations. 

a.  Bottom®,  @®Q  ©,  © 

yW.0.+*)  =  0 . (48) 


ci  •' 


d2  dd+*  d* 

~ 2  V+*  T*  V  *  d  \ 
C1  1 
b.  Left  Side®®® l© 


d$  .  .  .  .(49) 
=  0  .  .  .(SO) 


=  0 


(51) 


2  ee,  e. 

x  +  x  +— i.  x  -  JL,  -  0 

<M>  _2«W>c24<ei> 

C1  C1 

.  e*  e* 

yw  +  ~y«+"iyr~ 

ci  ci 


y .  =  0 

e  r<t> 


.(52) 

.(53) 


Right  SideQQQ© 


r((J,+,Nj)  =  constant . (54) 

for  x  the  same  as  Eq.  52 
for  y  the  same  as  Eq.  53 


d.  Upstream  EntranceQ^Q),  @.00® 


x(o,+.+  *)  =  o . (55) 

y(0.+,+*)  =  ^  H . (56) 

=  j-jj-  H . (57) 


e.  Downstream  Exit®®  @.0®  © 

=  constant  . (58) 

for  y  the  same  as  Eq.  56 
for  z  the  same  as  Eq.  57 

f.  Free  SurfaceQ )®  @00.® 


(g2+  i2> 


(eg  +  ih) 


|x 


(e2+h2) 


+  c 


2 

1 


2agN2D2(H- 


y) 


o  .  .  .(59) 


x(tf>.  = 


V  NjD  (2agMH-y)  V  w  V 


(60) 


•l 


dx  dy  Oz 

C1  dh  +  §+*  64* 


.  .(61) 

g.  Le£t  Side  of  Strut  @.  0 

s(4,+,NS)  =  f(<h)  specified  by  input  .  .(6?) 


2  eej  e. 

x.  .+  x  +  — *  x - -  -  0  .  .(63) 

<M>  2  +  e 

C1  C1 


y  +-£iy  +  _£y  _  *  =  o  .  .(64) 

_  2  2  y*  e  y<t> 

1  C1 

h.  Rj^ht  Side  of  Strut  @.  (g),  @,  @  @ 
for  z,  x,  and  y  the  same  as  Eqs.  62, 63,  and  64 


A  number  of  the  boundary  conditions  just  given  are 
immediately  obvious.  The  equations  for  other  boundary 
conditions  result  only  after  some  algebraic  manipulation. 
For  these  latter  conditions  an  explanation  and  the 
derivation  of  the  given  equation  are  contained  below. 


The  condition  for  x  along  the  channel  bottom  has 
been  obtained  from  Eq.  20  by  noting  that  since  y  =  0 
(constant)  along  this  boundary  y<j,*  =  0.  Therefore, 


cix*  =  y+v . (49a) 

Integrating  Eq.  49a  with  <J»  and  4>*  both  held 
constant  gives  Eq.  49  in  which  the  subscripts  to  the 
integral  sign  denotes  that  ifi  and  are  to  be  held 
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constant.  Since  is  constant  along  the  bottom,  the 
integration  of  Eq.  49  becomes  a  numerical  line  integration 
along  the  bottom  \p *  =  constant  grid  lines  with  $  being 
incremented.  The  implementation  of  this  line  integration 
involves  the  two  step  procedure  of  first  evaluating  tire  two 
partial  derivatives  with  the  integral,  and  secondly  carrying 
out  the  numerical  integration.  The  evaluation  of  y.  has 
been  based  on  up  through  third  order  forward  differences 
and  is  given  by 


S+I-.J3  A  4*  [3yi2k  *  6  yilk  *  2  Vi3k  +  3  yi4k]  *^6S^ 

The  equation  of  z^*  has  been  based  on  third  order 
forward,  fourth  order  central,  or  third  order  backward 
differences  respectively  depending  upon  whether  *J>*  =  0, 
4>*  =  &  ip*,  whether  lies  within  the  central  portion  of 
the  bottom,  or  whether  i|»*  =  N, ,  or  i|»*  =  M,  -  AiJ/*.  For 
tji*  =  0  (i.e.  k  =  1)  the  equation  is 

dz  .  if,  !i  3  .1  1 

§+♦  I k_ ! *  L  *  6  *11,1*2  *it, 3  3  *il,4j 


The  equation  for  =  Aiji*  =  1  is 


jz  |  1  r  l  l  \  I 

d+*|k=2  A4»*  [*il,3  *  3*»l,t*  2  *il,2  *  6  *iMj 


For  the  central  portion  (i.e.  k  =  3, 4, ...,  N-?) 


ct  |  i  r  2 ,  ,  i  , 

5+*}“  A+*  I  3  *ilk+ 1  ”  *ilk-l  *  12 


ilk+  2  *  *ilk- 


Equations  similar  to  Eqs.  67  and  68  but  based  on 
backward  differences  apply  along  the  lines  4'*  =  N,  (k  = 
N)  and  =  -  At|>*  (k  =  N-l ). 

The  numerical  integration,  for  other  than  the  first 
interval  <f>  =  0  to  $  =  =  1  or  the  final  interval  $  =  L, 

■  A$  to  ♦-  L ,  has  been  based  on  Bessel's  interpolation 
formula  for  a  third  degree  polynomial  passing  through  the 
values  of  four  consecutive  values  of  the  arguments  given 
by  the  product  of  the  two  derivatives  and  z^* .  This 
integration  formula  is 


Across  the  -first  and  final  intervals  this  integration  has 
been  based  on  the  trapezoidal  formula, 


A  <t> 

=  xLv 


+  y*V|. 


in  which  i  =  1  or  i  =  L-l. 

The  other  boundary  conditions,  Eqs.  60  and  61. 
which  contain  integrals  as  in  Eq.  49  have  been  handled  in 
the  same  general  manner.  Individual  details  differ  in  each 
case,  but  the  derivatives  involved  in  the  particular  equa¬ 
tion  have  been  approximated  by  third  order  forward  or 
backward  difference  and  fourth  order  central  difference  in 
the  interior  wherever  possible.  The  integration  has  been 
based  on  a  polynomial  passing  through  values  equal  to 
four  derivatives  or  products,  etc.,  thereof,  except  over  the 
first  and  last  intervals  of  the  line  integration  which  has 
been  based  on  the  trapezoidal  rule.  If  the  final  value  of 
the  integrated  variable  is  known  upon  reaching  the 
opposite  boundary  (as  with  z  from  Eq.  61  on  the  free 
surface)  any  error  in  not  closing  on  this  correct  value  is 
proportioned  according  to  the  distance  from  the  begin¬ 
ning  point.  In  the  case  of  x  on  the  bottom  and  the  free 
surface  (Eqs.  49  and  60)  the  average  of  all  final  values  is 
obtained  and  then  the  individual  differences  from  this 
average  are  distributed  according  to  the  magnitude  of  x. 

The  boundary  condition  for  z  along  the  bottom,  Eq. 
50,  has  been  obtained  from  Eq.  32  by  noting  that  f  =  y^* 
=  0  and  consequently  Eq.  32  reduces  to  Eq.  50. 

The  finite  difference  operator  for  Eq.  50  is: 

-(1  +  03,*i-l,lk+2(1  +  0l,*Uk-(1-°3,*i+Mk 

*  (Va2,*ilk-l+(V°2,*ilk+l  .  .  -  -(71) 

in  which 

ai=  V2'Ct=  (yi2k*yilk,2/cl 

°2  =  2  y+y4-^/Cl  =  4  (yi2k*yilk,(yi2k+ryilk+l 
"yi2k-l  +  yilk-l,/cl 

°3  *  2  VV  4  (yi+l,  2k*  yi+l,2k*  yi-l,2k 
+  yi-1.2k,/(yi2k-yitk) 

"i.ic  boundary  condition  Eqs.  52  and  53  which 
apply  for  x  and  y  respectively'  on  the  two  sides  of  the 
channel  are  obtained  irom  Eqs.  23  and  24  by  noting  that 
z  is  constant  in  the  planes  of  the  channel  sides  and 


consequently  S  =  z\b  =  0  and  c  =  z<j>  =  0.  The  finite 
difference  operators  for  Eqs.  52  and  S3  are  respectively 


{1  +  a3)Vljk1+2(1  +  VXijk1-(1-a3)Xi  +  ljk1 


|1  +  Vl.«k1t«U^-l,-,3W 

=  (VVVlk,  +  (Va2>W  .  -(73) 

in  which 

e2  ,2.  ? 

a.  =  — r  =  (z.  -  z  )  /c, 

1  _  2  ljk.  ljk,  1 


°2  2  "  8  ^ijkj'^ijk^^ij+lkj'^ij-lkj 

+  z. .  ..  )/c  2 
>J-lk2  1 

*<J  t 

°3  -  2e  -  4  ^*i+ljkj  *i+  ljk2  Zi-ljkj +  'i-ljk^ 

<Z‘Jk, "  *Uk2) 

and  the  subscripts  k,  =  1,  k2  =  2  on  the  left  side  and  k,  = 
N,  k2  -  N-l  on  the  right  side  of  the  channel. 

At  the  upstream  entrance  the  boundary  conditions 
have  been  obtained  under  the  assumption  that  the  flow  is 
uniform.  The  value  of  x  is  therefore  constant  and  assigned 
the  value  of  zero,  and  the  difference  in  z’s  between 
consecutive  <l>*  constant  grid  lines  on  the  face  equals  the 
difference  in  y’s  between  consecutive  <p  constant  grid 
lines.  The  width  of  channel  is  therefore  specified  in 
relationship  to  the  depth  of  flow  by  the  equation 


solved  using  only  1/2  of  the  space.  The  full 

problem  was  used  primarily  to  examine  the  performance 
of  the  solution  technique  in  duplicating  the  symmetry 
which  must  exist.  For  a  nonsymmetry  problem,  the 
specification  of  the  shape  of  the  strut  could  be  given  by 
changes  of  z  from  that  at  the  stagnation  equal 
constant  line  at  the  front  of  the  strut. 

The  boundary  conditions,  Eqs.  63  and  64,  for  x  and 
y  on  the  sides  of  the  strut  are  obtained  by  noting  that  z  is 
only  a  function  of  <j>  but  constant  with  variation  of 
Consequently  g  =  z^  =0,  and  Eqs.  and  24  reduce  to 
Eqs.  63  and  64,  just  as  they  reduced  to  Eqs.  52  and  53  for 
the  boundary  conditions  on  the  channel  sides.  Since  Eqs. 
63  and  64  are  identical  to  Eqs.  52  and  53,  the  finite 
difference  operators  for  x  and  y  along  the  strut  are 
identical  to  Eqs.  72  and  73  with  k,  and  k2  equal  to  NS 
and  NS-1  respectively  along  the  left  side  of  the  strut  and 
k ,  =  NS,  k2  =  NS+1  on  the  right  side  of  the  strut. 

Free  surface  boundary  conditions 

The  implementation  of  the  previously  discussed 
boundary  conditions  has  been  relatively  straight  forward. 
This  has  not  been  the  case  with  the  conditions  on  the  free 
surface,  and  consequently  the  free  surface  conditions  and 
the  difficulties  encountered  in  implementing  them  will  be 
discussed  in  this  section. 

The  free  surface  is  at  atmospheric  pressure  and 
consequently  the  sum  of  the  velocity  head  and  the 
magnitude  of  y  at  any  point  on  the  surface  must  equal  a 
constant  Hg ,  i.e. 


w  =  m-;h 


in  which  ag  is  the  acceleration  of  gravity.  The  magnitude 
of  the  velocity  squared  can  be  expressed  in  terms  of  the 
inverse  variables  by  substituting  each  of  the  first  equations 
in  Eqs.  6,  7,  and  8  into  the  middle  portions  of  Eqs.  1.2, 
and  3  respectively  to  give. 


The  assumption  of  uniform  flow  is  also  made  at  the 
downstream  final  face.  Therefore  the  conditions  are  the 
same  as  at  the  upstream  face  except  that  x  becomes  the 
constant  equal  to  the  average  of  the  x’s  obtained  from  the 
numerical  integrations  of  Eqs.  49  and  60  on  the  bottom 
and  free  surface  respectively. 


r2  =  u2+v2+w2=J2^2+y,2+z/l  •  •  •  -(75) 


But  from  the  definition  of  J  below  Eq.  6  it  can  be  shown 
that  J  =  VJ,  therefore  Eq.  75  reduces  to. 


The  shape  of  the  strut  has  been  specified  by  giving 
the  z  coordinates  of  both  sides  along  a  specified 
constant  line.  In  the  example  solution  given  later  this  ip* 
constant  line  has  been  selected  midway  between  the  two 
channel  sides,  and  the  z’s  have  been  specified  for  each 
potential  surface  intersecting  with  both  sides  of  the  strut 
to  result  in  symmetry  about  the  center  plane  of  the 
channel.  This  simplified  problem  could  actually  have  been 


V  ♦  y  i  ♦  ^ 

Upon  combiring  Eqs.  76  and  74  and  transforming  the 
resulting  equation  by  Eq.  17  to  introduce  the  dimen¬ 
sionless  dependent  variabtc  4>  leads  to  the  following  basic 
free  surface  boundary  condition. 


p 

5- 

I 


i 

i 


P 


i 


i 


& 


i 


2  2  .  2 
X0  +  y0  +  *♦ 


2igNj2D2(H  -  y) 


(77) 


A  number  of  alternative  schemes  for  use  of  Eq.  77 
in  the  finite  difference  solution  for  establishing  the  values 
of  x.  y  and  z  on  the  free  surface  have  been  investigated 
by  implementing  them  in  computer  algorithm  and  then 
examining  their  performance.  All  of  the  schemes  which 
have  been  investigated  are  associated  with  a  number  of 
numerical  difficulties  and  it  appears  that  successful 
solutions  can  be  obtained  only  if  the  computational 
procedures  anticipate  and  avoid  these  potential  difficul¬ 
ties. 


The  first  scheme  implemented  developed  the  finite 
difference  operato>  for  y  on  the  free  surface  by  replacing 
the  first  term  in  Eq.  77  by  its  equivalent  from  Eq.  20  and 
assumed  the  third  term  might  be  considered  the  known  c 
=  z^.  The  derivatives  in  the  resulting  equation  were  then 
approximated  by  second  order  central  differences.  The 
value  yjmt  l  k  at  the  nonexistent  grid  point  in  this 
difference  equation  was  eliminated  by  combining  the 
equation  with  the  operator  Eq.  24  for  y  at  interior  grid 
points  to  give  the  operator. 


£  -  -<1  +  a3>yi-lmk+2(1  +  <Xt,yimk-‘1-V.- 
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*  2aiyim-lk  ~  |*  + 


3  i+  lmk 
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'r^V-w 


~£*i+lmk  *i-lmk^  "£yH-bnk"yi-lmk)2 


-°4  =  ° 


(78) 


in  which  the  o’s  are  as  defined  in  Table  2  for  Eq.  24. 

Note  that,  Eq.  78  cannot  be  solved  explicitly  for 
yhnk  even  under  the  assumption  that  the  derivatives  given 
by  the  single  letters  are  known  as  is  the  case  with  the 
other  operator  in  Table  2.  Consequently  in  applying  Eq. 
78  at  consecutive  grid  points  along  the  line  defined  by 
incrementing  i  in  the  LSOR  method  a  system  of  nonlinear 
algebraic  equations  results.  This  nonlinear  system  has  been 
solved  by  the  generalized  Newton  iterative  method, 


7(q+n  =  y(q)  D-lT(q)  . (79) 

in  which  "jT is  the  vector  of  unknowns  yin)k,  i  -  1, 2, ...  L* 
the  superscript  q  denotes  iteration  number,  the  vector  f 
has  as  its  elements  the  values  obtained  by  applying  Eq.  78 
at  the  individual  grid  points,  and  D  is  the  commonly 
defined  Jacobian  derivative  determinant. 
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Since  the  determinant  D  is  tridiagonal  for  the_system  of 
equations  obtained  from  Eqv78^and  since  D'lf  represents 
the  solution  to  a  system  Dx  =  f,  the  same  algorithm  used 
for  solving  the  other  operators  in  the  LSOR  method  can 
be  utilized.  The  difference  is  that  the  algorithm  must  be 
repeated  for  each  iteration  in  obtaining  the  solution 
across  any  line. 

In  studying  the  performance  of  this  approach  it 
soon  became  evident  that  it  required  a  very  good 
initialization  of  values  if  it  were  to  give  correct  values  for 
y  on  the  free  surface.  Not  only  did  it  require  a  close 
initialization  of  surface  values,  but  the  interior  values  had 
to  be  well  settled  by  their  operators  before  the  free 
surface  values  were  adjusted  Much  could  be  written  about 
why  this  is  necessary,  but  perhaps  the  best  means  for 
illustrating  why  good  initialization  is  critical  to  prevent 
the  free  surface  values  from  straying  too  far  from  their 
correct  values  during  the  process  of  obtaining  a  finite 
difference  solutic  r  is  by  graphs.  Figures  2A  through  2S 
show  the  variation  of  the  function  f  defined  by  Eq.  78 
over  a  relatively  small  range  of  yimk  In  each  case  H  was 
equal  to  10.5  and  the  remaining  values  in  Eq.  78  (with  the 
exception  of  perhaps  one)  were  assigned  the  correct  values 
for  uniform  flow  with  a  depth  y  =  ID.  Figure  2A  shows  f 
vs.  yimk  with  all  other  values  correct,  and  as  expected  f  is 
zero  when  yimk  =  10.  However,  an  additional  zero  of  f 
exists  for  a  slightly  larger  value  of  yjmk  also.  Obvious 
numerical  difficulties  would  result  if  the  incorrect  second 
zero  for  f  were  selected  during  the  iterative  solution 
process. 

More  critical,  however,  is  the  behavior  of  the 
function  with  small  changes  of  the  other  variables  in  the 
equation.  The  remaining  curves  in  Figs.  2  show  ihe 
variation  of  f  with  y^  with  one  of  the  other  dependent 
variables  x,  y,  or  z  at  an  adjacent  grid  point  set  to  a  value 
±  0.15  from  its  correct  values.  As  the  value  yj+imk 
decreases  slightly  from  its  correct  value  (see  Fig.  2B),  it 
causes  the  first  zero  of  f  to  correspond  to  a  value  of  y!mk 
less  than  10  as  one  would  expect  and  the  second  zero  to 
become  further  removed  from  the  first.  With  an  increase 
in  yj4  lmk,  however,  the  two  zeros  became  closer  together, 
and  for  the  increase  of  0.15  used  in  Fig.  2C  to  almost 
coincide.  Cleatly  a  larger  error  would  cause  no  real  zero  to 
exist  for  the  function  at  least  in  the  vicinity  of  the  value 
yimk  needed  f°r  a  correct  solution.  When  this  s:*uation 
occurs,  the  Newton  method  would  project  off  far  from 
the  value  being  sought.  This  grossly  erroneous  value  in 
turn  affects  adjacent  functions  with  the  result  that  the 
Newton  method  never  returns  during  subsequent  iter¬ 
ations  with  reasonable  values  ■'f  y^v 
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which  must  be  satisfied.  In  implementing  Eq.  59  for  y  it 
was  decided  that  its  finite  difference  equivalent  should  be 
satisfied  on  a  point  by  point  basis,  rather  than  over  an 
entire  line  at  a  time,  because  of  the  difficulties  associated 
with  coping  with  vanishing  zeros  of  the  function. 

The  boundary  condition  Eq.  59  for  y  on  the  free 
surface  is  obtained  by  substituting  for  m  and  from 
Eqs.  20  and  22  into  Eq.  77.  The  finite  difference  operator 
for  y  on  the  free  surface  has  been  obtained  by  replacing 
the  derivatives  in  Eq.  59  by  second  order  or  higher  order 
differences.  For  those  grid  points  with  k  =  2,  k  =  N,  or  i  = 
2,  the  operator  will  be  different  than  at  the  interior  point 
on  this  boundary  because  of  the  need  to  use  non- 
symmetric  differences  to  approximate  yp*  .  For  interior 
free  surface  boundary  grid  points  the  operator  is 


Fig.  2.  Continued. 


Decreasing  or  increasing  the  other  values  of  x.  y,  or 
z  at  surrounding  grid  points,  or  x  and  z  at  that  free  surface 
grid  point,  cause  similar  behavior  in  f  as  can  be  noted  by 
the  remaining  curves  in  Figs.  2.  For  the  small  errors  of  i 
0.1 5  used  in  plotting  Figs.  2  the  function  f  has  lost  its  real 
zeros  in  12  of  the  19  cases  shown. 

An  operator  for  x  developed  in  much  the  same 
manner  as  Eq.  78,  except  that  the  second  term  in  Eq.  77 
was  replaced  by  its  equivalent  from  Eq.  21,  and  the  result 
combined  with  Eq.  24,  was  also  studied.  It  w  .  hoped  that 
since  the  yimk  which  appears  in  Eq.  78  and  also  appears  in 
the  finite  difference  operator  for  x,  but  in  the  x  operator 
it  is  considered  known,  that  the  implementation  of  the 
operator  for  x  would  be  associated  with  fewer  problems. 
However,  as  with  the  y  operator  convergence  was  ob¬ 
tained  by  the  Newton-LSOR  methods  combined,  only  if 
the  initialization  was  very  good. 

Similar  problems  exist  for  an  operator  for  z  ob¬ 
tained  by  combining  the  finite  difference  expression  from 
forms  of  Eqs.  77  with  the  regular  z  operator.  It  has 
therefore  been  concluded  that  the  difficulties  discussed 
above  are  an  inherent  part  of  the  free  surface  boundary 
condition  Eq.  77  which  cannot  be  eliminated  by  some 
manipulation  of  the  available  equations.  At  least  no  such 
manipulation  is  obvious.  It  was  decided,  therefore,  that 
programming  techniques  would  have  to  be  adopted  that 
would  prevent  the  difficulties  from  occurring  as  far  as  this 
was  practical,  and  that  means  would  have  to  be  included 
in  the  program  to  handle  the  difficulties  when  they  did 
occur. 

The  approach  for  handling  the  free  surface 
boundary  conditions  for  x  and  z  was,  therefore,  based  on 
the  integral  Eqs.  60  and  61,  and  Eq.  59  was  used  for  y  on 
the  free  surface  to  incorporate  the  basic  condition  Eq.  77 
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For  grid  points  along  the  grid  line  with  k  =  2,  the  quantity 
within  the  first  square  brackets  has  been  replaced  by  the 
following  third  order  difference  approximation 


V 1  *  yim3  -  I  yiml  -  I  yim2  *  { yin>4  -  -  (81) 

•k=2 

and  for  the  grid  line  i  =  2.  Similar  changes  are  necessary  in 
the  second  from  the  final  term  in  Eq.  80. 

Before  describing  how  Eq.  80  has  been  solved  to 
overcome  the  difficulties  associated  with  the  vanishing  of 
the  real  zeros  and  dual  zeros  of  the  previous  operator  f  of 
Eq.  78,  it  is  well  to  point  out  some  characteristics  of  Eq. 
80.  First  note  that  Eq.  80  is  also  implicit  in  yjmk  and 
therefore  the  solution  to  Eq.  80  must  be  obtained  by  an 
iterative  method.  Plots  showing  the  variation  of  f  of  Eq. 
80  with  yi(nk  for  small  errors  in  some  of  the  surrounding 
points  are  given  in  Fig.  3.  These  graphs  indicate  the 
following: 
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1 .  The  zero  being  sought  is  always  associated  with 
the  smaller  of  the  two  possible  values  of  ybnk  which  lie 
within  the  vicinity  of  the  solution  being  sought.  (While 
the  smaller  yitnk  is  sought  for  the  formulation  given  earlier 
in  this  report,  it  does  not  indicate  this  will  always  be  so.  It 
was  discovered  that  if  both  the  dependent  and  indepen¬ 
dent  variables  are  made  dimensionless,  then  the  larger  of 
these  two  possible  roots  is  the  correct  one.) 

2.  The  function  has  a  positive  derivative  with 
respect  to  y^  at  the  point  of  the  correct  zero  if  this  zero 
exists. 

3.  The  function  has  a  negative  derivative  with 
respect  to  y^  at  the  point  of  the  incorrect  zero. 

4.  If  real  zeros  of  the  function  vanish,  then  yimk 
associated  with  the  minimum  absolute  value  of  the 
function  (or  the  minimum  of  the  function  squared)  is  not 
too  different  from  the  value  of  y^  which  is  being 
sought,  and  which  will  again  come  into  existence  when 
the  errors  in  the  variables  at  the  surrounding  grid  points 
are  eliminated. 


Based  on  these  observations,  the  solution  to  Eq.  80 
at  each  grid  point  on  the  free  surface  proceeds  at  first  by 
computing  both  f  and  its  derivative  based  on  the  current 
value  of  yimkand  also  for  a  value  yhnk+  Ay,  in  which  A  y 
is  small,  i.e.  A  y  =  .005.  The  derivative  of  f  is  given  by: 
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(It  is  necessary  to  modify  Eq.  82  slightly  for  grid  points 
adjacent  to  the  boundaries,  i.e.  for  k*2,  k  =  N-2,  and  i  = 
2.) 


Should  the  derivative  associated  with  yiwk  be 
negative  (obviously  the  other  derivative  associated  with 
yj,„k  +  Ay  will  also  be  negative),  then  ytak  is  decreased  by 
a  small  amount  (in  the  order  of  .005)  and  the  procedure  is 
begun  again.  If  upon  repeating  this  process  of  decreasing 
yimk  *1*  fits*  derivative  becomes  positive  and  the  second 
remains  negative  and  both  of  the  functions  are  still 
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negative  then  no  real  root  exists  (or  if  it  did  exist  two 
zeros  would  lie  between  ytak  and  yimk  +  .005)  and  the 
solution  proceeds  by  calling  on  an  algorithm  which  uses 
the  Fibonacci  search  technique  (Wilde  and  Beighter, 
1967)  to  find  the  yimk  associated  the  minimum  value  of 
f 2  within  the  interval  yimk  to  yimk+  .005. 


Should,  on  the  other  hand,  the  first  computation 
for  f  be  positive  and  the  derivative  3f/  3y  also  be  positive 
then  yimk  is  decreased  in  value  until  f  associated  with  yimk 
becomes  negative  and  f  associated  with  yimk  +  .005  is  still 
positive.  As  soon  as  this  occurs  the  Newton  method 
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or  the  False -Position  method  is  used  to  find  the  zero  of  f. 


The  remaining  possibility  is  that  the  value  of  f 
associated  with  ybnk  is  negative  (its  derivative  positive), 
but  f  associated  with  yjmk  +  .005  is  also  negative.  Should 
this  be  the  case  then  yimk  is  increased  in  value  by  a  small 
amount  and  the  computation  repeated  until  one  of  the 
cases  above  has  occurred,  and  the  smaller  y^  giving  a 
zero  f  has  been  obtained  or  the  yimk  associated  with  the 
minimum  value  of  f2  has  been  found  by  the  Fibonacci 
search.  The  procedure  is  followed  through  for  each  point 
in  the  interior  of  the  free  surface  boundary.  Should  the 
procedure  fail  to  find  the  correct  zero  for  f  within  a 
specified  number  of  repetitions  in  increasing  or  decreasing 
yimk,a  message  is  printed  to  this  effect  and  the  average  of 
the  surrounding  four  y’s  in  the  plane  of  the  free  surface  is 
used.  Also  each  time  it  becomes  necessary  to  resort  to  the 
Fibonacci  search  to  minimize  the  absolute  value  of  f2 
because  a  real  zero  has  become  nonexistent  a  message  to 
this  effect  is  printed.  It  has  been  observed  that  after  the 
interior  values  are  well  settled  and  the  value  x,  y,  and  z  on 
the  free  surface  are  also  nearly  settled  that  the  minimiza¬ 
tion  procedure  is  called  upon  only  rarely  and  then  only 
near  the  stagnation  points  at  the  front  and  rear  of  the 
strut. 


The  boundary  condition  operator  Eq.  60  for  x  on 
the  free  surface  is  also  based  on  the  basic  condition,  Eq. 
77.  This  equation  is  obtained  by  solving  Eq.  77  for  xa 
and  then  integrating  the  resulting  equation  by  holding  $ 
and  41*  constant. 


The  boundary  condition,  Eq.  61,  for  z,  however, 
comes  from  solving  Eq.  20  for  z^,*  and  subsequently 
integrating  with  4*  and  i{>  held  constant. 
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Ffc.3.  V»i*tioo  of  the  finite  difference  function  jiren  by  Eq.  80  orer  a  *n*D  ruige  of  y^k with  no  error  (A)  or  a  small 
error  at  several  smoandng  pH  poiat*  (B-Q). 


Computer  logic  in  obtaining  solution 

The  general  logic  followed  in  writing  the  computer 
program  used  in  obtaining  the  solution  is  illustrated  in  the 
gross  flow  chart  contained  as  Fig.  4.  While  it  is  not 
necessary  to  follow  the  exact  pattern  of  obtaining 
tentative  solutions  given  in  Fig.  4,  it  is  desirable  to  provide 
control,  through  data  input,  so  that  tentative  solutions  on 
boundary  planes  are  not  necessarily  obtained  during  each 
cycle.  From  the  above  discussion  of  difficulties  associated 
with  obtaining  the  tentative  solution  for  y  on  the  free 
surface  boundary  it  is  clear  that  it  would  not  be  expedient 
to  call  on  the  subroutine  which  obtains  this  tentative 
solution  during  the  first  few  cycles,  particularly  if  a  rough 
initialization  is  used.  Also  not  adjusting  those  boundary 
values  which  are  obtained  by  integration  (see  Eqs.  49, 60, 
and  61)  until  after  the  interior  values  have  been  initially 
settled  contributes  to  more  rapid  convergence  to  the  final 
solution.  Consequently  specifications  are  included  as  part 
of  the  input  data,  which  determine  during  which  cycles 
each  of  the  subroutines  will  be  called  to  obtain  the 
tentative  solution  from  the  boundary  condition  operators. 
The  present  program  first  obtains  all  tentative  solutions 
for  x,  which  are  called  for  during  that  cycle;  then  it 
obtains  the  tentative  solution  for  y  before  obtaining  the 
tentative  solutions  for  z.  Tentative  solutions  on  interior 
planes  are  obtained  prior  to  obtaining  the  tentative 
solutions  for  that  variable  on  boundary  planes. 

The  logic  followed  in  obtaining  each  of  the  tentative 
solutions  varies  slightly  dependent  upon  which  variable  is 


involved  and  whether  the  plane  is  interior  or  a  boundary 
plane.  For  tentative  solutions  on  interior  planes  and  those 
boundary  planes  in  which  the  LSOR  method  is  used  to 
obtain  a  solution  to  the  boundary  condition  equation,  the 
flow  chart  given  as  Fig.  S  shows  the  essential  logic.  On 
those  boundary  planes  on  which  the  tentative  solutions 
are  obtained  by  numerical  line  integrations,  the  logic  is 
such  that  the  required  integration  is  completed  a  line  at  a 
time  beginning  with  a  line  on  or  adjacent  to  an  edge,  and 
proceeds  until  it  has  been  completed  for  the  last  line  on 
the  opposite  edge  of  the  plane. 

For  use,  particularly  during  earlier  cycles,  the 
program  provides  that  the  free  surface  values  of  y  may  be 
smoothed  if  desired.  This  smoothing  is  accomplished  by 
fitting  the  values  of  y  along  i  =  constant  lines  (i.e.  with  k  = 
1 ,  2, ...  N)  on  the  free  surface  by  least  squared  regression 
analysis  to  an  equation  of  the  general  form. 


yimk  =  b0+blk+b2k2+b3Cos(k-NS-1,t^  •  (84) 

in  which  the  b's  are  the  coefficients  obtained  from  the 
regression  analysis.  After  fitting  the  y’s  along  each  such 
line  to  an  equation  of  the  form  of  Eq.  84,  they  are 
adjusted  to  satisfy  the  equation  exactly.  The  adjustment 
also  includes  the  y  on  the  two  channel  walls  and  the  two 
sides  of  the  strut. 
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Fig.  4.  Flow  chart  of  logic  used  in  obtaining  the  solution. 
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F%.  5.  Flow  chart  of  fcfic  used  in  computer  program  subroutines  which  obtain  the  tentative  solutions  by  the 
LSOR-Method. 
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The  final  solution  consists  of  the  magnitudes  of  x, 
y,  and  z  at  all  grid  points  within  the  ifnfiip*  space  used  to 
solve  the  problem.  Consequently  the  coordinates  are  given 
for  each  intersection  of  the  potential  surfaces  with  all  of 
the  orthogonal  stream  surfaces  defined  by  holding  ip  and 
ip  *  equal  to  constants.  In  this  form  the  solution  is  ideally 
adapted  for  presentation  as  a  space  flownet.  Such  a  space 
flownet  is  constructed  by  simply  connecting  all  consecu¬ 
tive  points  defined  by  the  x,  y,  and  z  values  given  at  each 
grid  point  throughout  the  <j>^*  space  by  lines  in  an 
isometric  drawing  (or  other  graphical  projections  which 
show  depth  into  the  paper  as  well  as  the  shape  within  the 
pl:.ne  of  the  paper).  The  small  planes  defined  by  these 
lines  represent  the  sides  of  each  element  of  the  flownet. 
The  intersection  of  the  ip  and  ip*  constant  planes  define 
the  streamlines  of  the  flow.  The  velocity  is  inversely 
related  to  the  area  of  the  square  formed  by  the  tj ;  and  ip* 
equal  constant  lines  and  the  distances  between  consecu¬ 
tive  equipotential  surfaces  as  given  by  combining  Eqs.  10, 
11,  and  12  with  Eq.  76  in  various  ways.  (Equations  for 
the  velocity  and  its  direction  are  given  later.)  That  is  the 
velocity  is  greater  in  regions  in  which  the  volume  enclosed 
within  individual  cubes  (or  parallelepiped  elements  if  A<j> 
=  Ai p  =  Aip*),  of  the  flownet  is  smaller  than  in  those 
regions  in  which  this  volume  is  greater. 


While  a  complete  isometric  space  flownet  can 
readily  be  obtained  by  use  of  a  computer  driven  plotter, 
the  numerous  lines  resulting  therefrom  would  make 
visualization  of  the  complete  flow  difficult.  Alternatives 
are  to  plot  only  a  few  of  the  flownet  lines,  or  to  plot  only 
the  flownet  lines  in  key  planes.  Fig.  6  has  been  prepared 
by  using  this  latter  type  of  plot,  in  which  the  plane 
flownets  from  the  top,  rear,  and  right  side  are  given  in  an 
isometric  projection  of  the  problem. 


The  more  essential  specifications  used  in  obtaining 
the  solution,  whose  flownet  is  given  in  Fig.  6,  are  as 
follows:  (1)  The  depth  of  uniform  flow  upstream  from  the 
strut  equals  10  feet  (2)  The  number  of  (pip*  grid  planes 
equals  the  number  of  $ip  grid  planes  and  consequently 
the  width  between  channel  sides  is  also  10  feet.  (3)  The 
number  of  ipip*  planes  (increments  in  the  $  direction 
phis  one)  was  given  as  20,  resulting  in  a  length  from 
beginning  to  end  of  the  problem  equal  to  18.4  feet.  (4) 
The  strut  was  specified  0.6  feet  wide  at  its  widest  point 
and  it  began  on  the  7th  ip  if/*  plane  and  ended  on  the  14th 
ifMp*  plane  resulting  in  a  length  equal  to  6.4  feet.  (S)  The 
velocity  head  in  the  undistributed  uniform  flow  equals  0.5 
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feet  (H  -  10.5  feet),  resulting  in  an  upstream  velocity 
equal  to  5.675  fps.  In  solving  this  problem  2,420  finite 
difference  grid  points  were  used.  Since  three  unknowns 
must  be  solved  for  simultaneously,  three  times  this  many 
finite  difference  grid  points  or  7,260,  were  actually  used. 

The  solution  to  this  problem  was  obtained  in  a 
piecemeal  manner  as  the  separate  subroutines  were  de¬ 
bugged,  etc.  Consequently  it  is  not  possible  to  give  the 
exact  amount  of  computer  execution  time  required  for 
the  final  solution.  With  an  initialization  which  is  easily 
generated  in  a  computer  program,  and  using  the  number 
of  grid  points  used  for  this  problem,  a  reasonable  estimate 
of  the  execution  time  on  an  UNIVAC  1108  system  is  15 
minutes,  however. 

While  an  isometric  drawing  of  the  space  flownets 
helps  in  visualizing  the  complete  flow  process,  more 
detailed  information  regarding  special  features  of  the  flow 
can  be  obtained  by  examining  the  flow  in  separate  planes 
within  the  space.  The  solution  from  an  inverse  formula¬ 
tion  is  in  an  ideal  form  to  examine  the  flow  field  in 
separate  equipotential  planes,  i.e.  planes  defined  by  ip  and 
ip*  axes,  or  for  examining  the  behavior  of  the  flow  in 
separate  planes  defined  by  <p<p  or  <P  <P*  axes.  A  solution 
to  a  three-dimensional  problem  in  the  physical  space  (i.e. 
in  the  space  defined  by  the  x,  y,  z  cartesian  coordinates) 
would  be  well  adapted  for  examining  details  in  separate 
xy  planes  (i.e.  defined  by  z  equal  a  constant),  xz  planes  or 
yz  planes  but  would  require  interpolation  to  examine  the 
flow  field  in  equipotential  planes  for  instance.  On  the 
other  hand  the  results  from  the  inverse  solution  in  the 
<p  ip  ip*  space  require  interpolation  to  examine  or  display 
the  flow  in  separate  planes  of  the  physical  space.  Thus  for 
example,  if  one  wishes  to  examine  the  flow  field  in  an  xy 
plane  with  z  equal  to  a  given  constant,  it  would  be 
necessary  to  obtain  the  magnitudes  of  x  and  y  which 
define  the  intersection  of  the  plane  flownet  lines  by 
interpolation  of  the  x’s  and  y’s  on  the  two  adjacent 
inverse  planes  that  contain  z  values  which  bracket  the 
specified  constant  z.  Obviously  accomplishing  this  is  not 
difficult;  perhaps  even  less  difficult  than  plotting  a 
flownet  given  a  solution  of  the  potential  function  in  the 
physical  space.  However,  no  flownets  from  such  planes  \ 

within  the  physical  space  are  given  herein.  For  boundaries 
on  which  either  x,  y,  or  z  is  constant  such  as  the  sides,  or  \ 

beginning  and  end  of  the  channel  problem,  no  interpola-  X 

lion  is  necessary.  The  flownets  from  such  boundaries  are  J 

simultaneously  on  a  plane  in  which  <P,  ip  or  ip*  is  J 

constant  as  well  as  x,  y,  or  z  is  constant. 
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Isometric  drawing  of  space  flownet  of  a  strut  in  a  free  surface  channel  flow. 
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In  general,  however,  individual  plane  flownets  ob¬ 
tained  by  holding  41,  \p,  or  ip*  constant  are  more  useful  in 
visualizing  flow  patterns  than  those  obtained  by  holding 
x,  y,  or  z  constant.  In  understanding  plane  flownets 
obtained  by  holding  4> ,  ip,  or  ip*  constant,  one  needs  to 
interpret  the  results  as  projections  upon  a  specified  plane. 
Thus  for  example  if  one  wished  to  examine  the  flow 
pattern  along  the  sides  of  the  strut,  a  plane  flownet 
obtained  by  plotting  x  and  y  from  the  plane  with  iJj* 
equal  to  the  value  coincident  with  the  strut,  would  be  a 
projection  of  the  strut’s  flow  pattern  upon  a  vertical  plane 
parallel  to  the  sides  of  the  channel.  Such  a  flownet  is  given 
in  Fig.  7. 


in  which  c,  =  D/M,  is  as  previously  defined,  Vc  is  the 
undistributed  upstream  velocity,  and  A$  is  unity  for  the 
solution  given  herein  and  can  be  deleted  from  the 
equation. 


The  direction  angles  of  the  velocity  vector  at  any 
point  within  the  flow  space  can  be  obtained  by  first 
noting  that  the  first  equation  of  Eqs.  6,  7,  and  8  can  be 
written  respectively  as: 


v  *♦  -  « 


A  plane  flownet  obtained  by  plotting  the  y  and  z 
magnitudes  of  the  solution  from  a  |  equal  constant 
(equipotential  plane)  is  given  in  Fig.  8.  This  flownet 
represents  a  projection  of  the  flow  from  the  leading  edge 
of  the  strut  unto  a  vertical  plane  at  right  angles  to  the 
channel  sides. 


V  y*  =  -v  . 


V  z.  =  w 


The  plane  flownet  obtained  from  the  free  surface  ip 
equal  constant  plane  is  shown  in  Fig.  9,  in  which  the  x 
and  z  magnitudes  have  been  plotted.  This  flownet  can  be 
interpreted  as  a  projection  of  the  free  surface  flow  pattern 
into  a  horizontal  plane  (i.e.  a  plane  parallel  to  the  channel 
bottom). 


These  equations  are  obtained  by  noting  that  the  Jacobian 
determinant  J  equals  the  magnitude  of  the  velocity 
squared.  Since  u  =  V  cos  a,  v  =  V  cos  B,  and  w  =  V  cos  y 
(in  which  a,  B,  and  y  are  the  angles  of  the  direction 
cosines  for  the  velocity  vector),  Eqs.  86,  87,  and  88  give 


a  =  cos  (V  x^) 


The  magnitudes  of  x,  y,  and  z  at  each  intersection 
of  grid  planes  in  the  space,  which  constitute  the 

basic  solution  can  readily  be  used  to  obtain  other 
quantities  of  interest  about  the  flow.  The  local  velocity 
magnitude  can  be  computed  from  the  following  finite 
difference  equation  derived  from  Eq.  76. 


P  =  cos'  (Vy^) 


y  =  CO.  (Vz^) 


2A(tc,V 
1  o 


The  pressure  at  any  point  can  be  computed  from 
Eq.  85  and  the  Bernoulli  equation,  i.e. 


*  fxHljk‘Xi-ljk)  +  (yi+ljk'yi-ljk)  +  (zi+$k*zi-ljk)Z] 


Pijk  =  pg(H'yijk)  "p  2  .  .  .  .(92) 


Fig.  7.  Plane  flownet  from  the  plane  associated  with  k=6  which  coincides  with  the  stream  surface  of  the  strut 
obtained  by  projecting  the  nugnhodesof  x  and  y  onto  a  vertical  plane  parallel  to  the  channel  sides. 
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F&.  8.  Plane  flownet  from  the  ty*  plane  associated  with 
i-7  which  touches  the  leading  edge  of  the  strut  and 
which  results  from  projecting  the  magnitudes  of  y 
and  z  onto  a  vertical  plane  at  right  angles  to  the 
channel  sides. 
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Fig.  9.  Plane  flownet  from  the  plane  associated  with  j=l  1  which  coincides  with  the  free  surface  obtained  by 
projecting  the  magnitudes  of  x  and  z  onto  a  horizontal  plane  parallel  to  the  channel  bottom. 


34 


,  1*.  1/nij.i'j  k>i  t A-ii  W»*y.  '■»* l  K*.  £4»jk  «rHfcs.Wi rW<AJ?A»  W.VI  iUfAt , 


CONCLUSIONS  AND 


The  use  of  a  mathematical  formulation  which 
reverses  the  usual  role  of  variables  shows  promise  as  a 
valuable  tool  for  numerically  solving  certain  types  of 
three-dimensional  potential  fluid  flow  problems.  Like 
most  new  approaches,  however,  the  merits  of  the  methods 
need  to  be  further  evaluated  and  improved.  The  methods 
and  techniques  used  in  this  report  for  solving  the  inversely 
formulated  space  boundary  value  problem  represents  an 
initial  approach  which  is  workable,  but  which  will  no 
doubt  be  streamlined  and  improved  upon  with  time. 

The  interchange  of  the  conventional  dependent  and 
independent  roles  played  by  the  variables  in  a  three- 
dimensional  potential  fluid  flow  problem  results  in  ad¬ 
vantages  similar  to  those  which  occur  in  solving  two- 
dimensional  plane  and  axisymmeiric  potential  fluid  flow 
problems.  Perhaps  the  major  advantages  are:  (l)That  the 
region  of  the  space  boundary  value  problem  is  a  paral¬ 
lelepiped  with  planes  for  boundaries,  which  in  the 
physical  plane  may  be  irregular  and  of  unknown  position, 
such  as  free  surfaces  or  cavity  surfaces,  and  (2)  the  form 
of  the  solution  is  better  adapted  for  graphical  presentation 
and  for  computing  various  items  of  interest  about  the 
flow.  These  advantages  occur  at  the  expense  of  more 
complex  simultaneous  partial  differential  equations. 

In  order  for  the  inverse  solution  method  to  be 
readily  adaptable  and  used  practically  for  solving  a  variety 
of  problems  involving  free  surfaces  and  cavities,  alternate 
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and  better  schemes  or  methods  are  needed  for  handling 
boundary  conditions  resulting  from  a  constant  pressure 
free  surface  under  the  influence  of  gravity.  The  approach 
used  herein  is  associated  with  a  number  of  difficulties 
which  no  doubt  will  become  progressively  harder  to  cope 
with  as  the  complexity  of  the  problem  increases.  Con¬ 
sequently,  a  problem  with  a  three-dimensional  cavity  and 
free  surfaces  would  represent  a  difficult  undertaking 
without  better  methods  for  handling  such  free  surface 
boundary  conditions.  With  such  improved  methods,  the 
inverse  formulation  should,  in  fact,  provide  a  practical 
numerical  solution  procedure  for  solving  three- 
dimensional,  steady-state,  free  surface,  and  cavity  poten¬ 
tial  fluid  flow  problems. 


Even  if  more  satisfactory  methods  for  handling  free 
surface  boundary  conditions  are  not  developed,  the 
methods  still  represent  a  valuable  tool  for  solving  three- 
dimensional  problems  without  free  surfaces,  particularly  if 
the  problem  is  a  design  problem  instead  of  an  analysis 
type  problem.  In  a  design  problem  shapes  of  confining 
structures  are  sought  which  give  some  desired  flow 
characteristics.  The  inverse  formulation  is  particularly  well 
adapted  for  such  problems  in  which  the  shape  of  a 
boundary,  which  is  a  stream  surface,  is  part  of  the 
solution,  resulting  from  a  specification  of  fluid  behavior, 
but  less  well  adapted  if  non-plane  confining  surfaces  have 
specified  shapes  as  in  analysis  type  problems. 
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